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Context. With the Herschel Space Observatory, lines of simple molecules (C + , O, and high-/ lines of CO, / up > 14) have been 
observed in the atmosphere of protoplanetary disks. When combined with ground-based data on [CI], all principle forms of carbon 
can be studied. These data allow us to test model predictions for the main carbon-bearing species and verify the presence of a warm 
surface layer. The absence of neutral carbon [CI], which is predicted by models to be strong, can then be interpreted together with 
ionized carbon [C II] and carbon monoxide. 
Aims. We study the gas temperature, excitation, and chemical abundance of the simple carbon-bearing species C, C + , and CO, as 
well as O by the method of chemical-physical modeling. Using the models, we explore the sensitivity of the lines to the entering 
parameters and constrain the region from which the line radiation emerges. 

Methods. Numerical models of the radiative transfer in the lines and dust are used together with a chemical network simulation and 
a calculation of the gas energetics to obtain the gas temperature. We present our new model, which is based on our previous models 
Q ■ but includes several improvements that we report in detail, together with the results of benchmark tests. 

Results. A model of the disk around the Herbig Be star HD 100546 is able to reproduce the CO ladder together with the atomic fine- 
structure lines of [O I] and either [C I] or [C II]. We find that the high-/ lines of CO can only be reproduced by a warm atmosphere 
with r gas » 7d ust . The low-/ lines of CO, observable from the ground, are dominated by the outer disk with a radius of several 100 
AU, while the high-/ CO observable with Herschel-PACS are dominated from regions within some tens of AU. The spectral profiles 
of high-/ lines of CO are predicted to be broader than those of the low-/ lines. We study the effect of several parameters including 
the size of the disk, the gas mass of the disk, the PAH abundance and distribution, and the amount of carbon in the gas phase. 
Conclusions. The main conclusions of our work are (i) only a warm atmosphere with T m » Ti ulit can reproduce the CO ladder. 
( ii) The CO ladder together with [O I] and the upper limit to [C I] can be reproduced by models with a high gas/dust ratio and a low 
abundance of volatile carbon. These models however produce too small amounts of [C II]. Models with a low gas/dust ratio and more 
volatile carbon also reproduce CO and [O I], are in closer agreement with observations of [C II], but overproduce [C I]. Owing to the 
uncertain origin of the [C II] emission, we prefer the high gas/dust ratio models, indicating a low abundance of volatile carbon. 
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1 . Introduction violet (UV) radiation and X-rays. Through heating of the disk's 

upper atmosphere, a d isk-wind is driven, which may eventually 

£ ; At an important phase of a young low-mass star's life, the na- difjperse ±e difjk ( Alexander et aL 2006; Mano et al . 20 08; 

. tal envelope has dispersed, but the star is still surrounded by an |o^&Hollenbach 2009). The most abundant species tracing 

■ optically thick protoplanetary disk, which is thought to be the the ionizing radiation and warm tem perature of a few 100 K are 

location of planet formation (see | Dullemond & Monniei| | 2010| (ionized) atomic carbon (C> c+)> atomic oxygen (Q)> and high . y 

and | Armitage| | 201fj| for recent reviews). Our own solar system Hnes of CQ (with j > 14) A from neutral carbon> which 

contains both gaseous giants and rocky planets. It is thus evident . g accessible to ground -based telescopes, these species can now 

that both the dust and the gas content of the disk need to be stud- be observed using the PACS instmmen t (Poditsch et al. 2010) 

led in order to understand planet formation. The feedback of the onboard the Herschel Space observatory (Pilbratt et al. 20k3T 
forming star onto the disk can help decide the disk's fate. The 

star irradiates the disk, heating both the gas and dust and chang- The main forms of carbon ( C > C+ ' and CO ) can thus be ob " 

ing the chemical composition of the gas by high-energy ultra- served for the first tlme and through high-/ CO lines also in 

warm regions. This allows us to probe the carbon chemistry 

Send offprint requests to: Simon Bruderer, and the physical structure of planet- and comet-forming zones, 

e-mail: simonbruderer@gmail.com which were previously inaccessible to us. The carbon budget 

* Herschel is an ESA space observatory with science instruments is important for planet formation, since the terrestrial planets 

provided by European-led Principal Investigator consortia and with im- in the inner Solar System are known to have a deficit in car- 

portant participation from NASA. bon of three orders of magnitude relative to the solar abundance 
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(Alle ereet all 1200 lh . Indications of a carbon deficit are also 
found towards other planetary systems (lJural2008h . In the diffuse 
ISM, about 40-60 % of the carbon is locked into amorphous car - 
bon grains dSavage & Sembachlll996l: ISofia et al.ll200l l201lh . 
To be consistent with the Earth's composition, these grains need 
to have been destroyed before planet formation. In comets, how- 
ever, carbon grains have been detec ted, suggesting t hat grain de- 
struction mechanisms play no role dLee et al1l2010h . 

What region of the disk does the far-infrared (FIR) C, C + , 
CO and O lines trace? How does the partitioning between 
volatile and refractory carbon vary from clouds to disks? The 
simultaneous measurement of the species tracing a significant 
fraction of the disk may provide answers to these questions. In 
addition, besides the hard to observe H2, these simple species 
form the most abundant constituents of the gas. Their lines are 
thus important tracers of the physical structure of the disk and 
allow us to finally test the predictions made by many exist- 
ing physical-chemical models of disks. Neutral carbon [CI], for 
example, is predicted to be st rong by models with a range of 
gas masses and dust settling (Jonkheid et al. 2007) but i s not 
detected i n observations of H P 100546 dPanic et al.l 2010) and 
CQ Tau (Chapillon et al. 2010). In addition to the submillimeter 
lines discussed in our work, neutral carbon also has forbidden 
lines in the near-infrared (at 8729, 98 26, 9852 A), which ar e 
predicted to have detectable strengths dErcolano et al.ll20 09b). 
To date, only a few detections towards pre-main-sequence stars 
have been reported (lLoope r et al. 2010). 

The physics and chemistry of gas in protoplanetary disks, 
which is irradiated by the central star, is similar to the in- 
terface of general molecular clouds, where nearby massive 
stars heat and photoionize the gas. Modeling the infrared 
(IR) response of gas irradiated by UV and X-rays has been 
developed in the context of these photo-dominated or disso- 
ciati on regions (PDRs) and X-ray dominated regions (XDRs) 
(e.g. iTielens & Hollenbach1ll985t ISternberg & Dalgarnol [19891 
iKaufman et al.l 1999; M eiierink & Spa ans 2005). These models 
simulate the infrared line emission by solving the coupled 
system of chemical evolution and thermal balance. The 
coupling is the result of the cooling rates for the thermal 
balance depending on the abundance of the coolants, which 
is determined by the chemical network that itself also de- 
pends on the temperature. These physical-che mical models 
were applie d in the cont e xt of protost ars by ISpaans et al.l 
(1 19951) and iBruderer et all d2009al l2010l) to the walls along 
outflows, irradiated by UV radiation arising from either 
accretion hot-spots or the protostar's photosphere. For pro- 
toplanetary di sks, physical-chemica l mod e ls have been 
applied by e.g. Kamp & van Zadelhoff d200lb. iGlassgold etafl 



(2004), 
d2004l). 



Nomura & Millar (2005), Aikawa & Nomura (2006), 
IJonkheid et all d2004l). IJonkheid et al.l (120061). IJonkheid et al.1 
d2007l). iGorti & Hollenbachl (120081)7 IWoitke et al I d2009a). 
IWoods & Willacy! d2009l) . and lErcolano et al.l d2009al) . 

The well-studied intermediate-mass pre-main-sequence star 
HD 100546 is the focus of our present study. The B9.5Vne 
type star of 2.2 M L hn , ~ 27 L , and T eS = 10500 K 
dvan den Ancker et al.l 119971) is at its distance of 103 ± 6 pc 
(Hipparcos) one of the closest Herbig Ae/Be sources. The 
star and disk have been extensively observed at all wave- 
lengths in lines and th e conti nuum. X-ray observations are re- 
ported in IStelzer et al.1 (2006). Both FUSE and IUE observed 
HD 100546 in the UV, detecting absorption in electronic tran- 
sitions of H 2 probing the warm gas with a temperature o f sev- 
eral 100 K dMartin-Zaidi et al.ll2008l) . lArdila et alJ d2007l) ana- 



lyzed scattered-light images in the optical using ACS/Hubble 
finding evidence of minimal grain sizes larger than those of 
typical ISM grains. The scattered light observations also reveal 
structures resembling spiral ar ms possibly caused by either the 
perturbations of a companion dOuillen et al1l2005h or a warped 
disk structure (Ouillen 2006). Features in the near-infrared 
(NIR) reveal a high cr ystalline silicate dust f raction in the in- 
ner disk surface layers dBouwman et alj|2003l) and an abundant 
amount of PAHs w hose emission is detected on extended scales 
dGeers et al 1 12007^ A wealth of continuum images in the near 
and mid-IR (e.g.lPantin et"al] |2000l: iGradv et al.ll200ll; iLiu et al.1 



2003; IGradv etaij f2005) have helped to constrain the structure 
of the inner disk using SEP fitting. They reveal a large inner 
hole of radius ~ 10 AU dMalfaitet all 119981; iBouwman etafl 
120031). which possibly has a Jupite r-sized planet orbiting in- 
side dAcke & van den Anckeri |2006|). Interferometric observa- 
tions (VLT- AMBER and MIDI. lLeinert et alj|2004l; iPetrov et afl 
2007) confirm t he presence of the i nner hole directly and set 
r hoIe ~ 13 AU dBenistv etafl 120101) . Some CO ro-vibrationa l 
emissi on at 4.7 was reported by Ivan der Plas et al. I d2009l) 
and lBrittain et all (|2009j). The latter paper inferred that the CO 
emission originates from radii of 13-100 AU, which is consistent 
with the inner hole seen in the dust continuum. They found a hot 
rotational temperature of ~ 1000 K and indications that both 
UV fluorescence and collisional excitation are important for the 
excitation. Hot gas in the disk surface layers is a lso detected by 
observing the H2 v = 1-0 S(l) line at 2.12 fj.m bv lCarmona et al.l 
d2011l) . 

The bulk of the disk mass is at cold temperatures, as traced 
by submillimete r lines and cont i nuum. Subm illimeter continuum 
was reported by Henning et al] dl994 H"998). who derived a dust 
mass of a few times 10~ 4 M Q . The gas up to ~ 100 K is traced 
by CO rotational lines up to J — 7 - 6, as reported bv lPanic et al.l 
(l2010l) . They derived a gas mass of > 10~ 3 M Q and concluded 
from the high 12 CO J = 6-5/ 12 CO 7 = 3-2 ratio that the disk at- 
mosphere needs to be warmer than in T Tauri stars, probably ow- 
ing to photoelectric heating of the gas by the strong UV field of 
the forming B9 st ar. Strong higher-./ fines up to CO J — 30 - 29 
were detected by ISturm et alll2010l using the PACS instrument 
onboard Herschel. The detected lines have upper level energies 
of up to 2500 K. 

In this work, we study the chemistry and excitation of sim- 
ple species (C, C + , CO, and O) in the disk around a pre-main- 
sequence Herbig Ae/Be star. Our approach is twofold. W e first 
want to interpret the PACS observatio ns of | Sturm et al] (2010) 
and the ground-based observations of iPanic et al.1 d2010l) taken 
towards HD 100546 and search for evidence of a warm disk at- 
mosphere from the CO ladder, which probes a wide range of 
temperatures. We however also wish to study the dependence of 
the lines on certain parameters such as the carbon abundance 
in the gas phase, the gas/dust ratio, and others. In Section [2] 
we briefly summarize the modeling approach and describe the 
observations. In Section [3] we first present the results for one 
particular model that agrees reasonably well with the observa- 
tions in Section [3] and then discuss the parameter dependences 
in Section|4]and|5] In Section|6] we summarize our results. 

2. Model and observations 



For this work, we use a new model based on IBruderer et al.l 
(2009b a, 2010), with several changes and improvements imple- 
mented. The section provides a short summary of the model. 
Details and results of benchmark tests are reported in Appendix 
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Our physical-chemical model starts with (i) a dust radiative 
transfer calculation for a given dust and gas distribution to obtain 
the local UV radiation field and the dust temperature. Using an 
initial guess of the gas temperature, ( ii) the chemical abundances 
are calculated, which serve as input for a molecular/atomic exci- 
tation calculation to obtain the cooling rates. Next, (Hi) the ther- 
mal balance is calculated to obtain an improved guess for the gas 
temperature. Steps (ii) and (Hi) are repeated until convergence is 
achieved. Raytracing and convolution to the telescope beam then 
yields line fluxes that can be observed with observations. These 
modeling steps are summarized in Figure lATl 

Differences from the study of lBruderer etal] (l2009bllaL[20Toh 
are that we solve the chemical network for the steady-state 
solution inst ead of interpolating fr om a pre-calculated grid of 
abundances dBruderer et alJ l2009b). This new approach allows 
us to include detailed photodissociation and ionization rates 
based on the local wavelength-dependent UV field and the ac- 
tual column densities for the self-shielding factors. The ther- 
mal balance is refined by obtaining atomic and molecular cool- 
ing rat es from the excitation calculation used in iBruderer et alJ 
(2010). T his (approxima te) method, which is similar to that of 
Poelman & Spaans (2005), allows us to calculate the molecular 
excitation quickly and with reasonable accuracy. We can now 
include the detailed velocity structure into the excitation calcu- 
lation and account for heating of the gas by the IR pumping of 
molecules, which is then followed by collisional deexcitation. 
The method couples cells at different radii of the disk, unlike the 
"vertical direction only" escape probability approach employed 
in most previous disk modeling. The disadvantage of our method 
is that we need to iterate over the entire chemical structure to 
derive the local gas temperature, which increases the calcula- 
tion time. However, the method allows us to calculate models 
with an arbitrary geometry including e.g. disks embedded in en- 
velopes. Since our code contains all geometrical properties in 
one module, it can be extended to three-dimensional structures in 
the future. The chemical network in the new code has also been 
extended by incorporating hydrogenation reactions on grain sur- 
faces, photodesorption of species frozen-out onto grains, and vi- 
brationally excited molecular hydrogen (Hp, which can be used 
to overcome activation barriers. 

In this study, we concentrate on the rotational ladder of CO 
and do not consider vibration-rotation lines detected in the M- 
band at 4.7 fim, since fluorescent excitation by UV photons plays 
a significant role in their population. We note that the excitation 
of the rotational ladder is determined by collisions and the rel- 
ative population of CO in the v = 1 state is small at 4 % for 
a thermalized population at 1000 K dBrittain et al.ll2009h . Thus, 
the results presented in this work would not change by includ- 
ing vibrationally excited levels of CO. For the same reasons, we 
also do not model the H2 2.12 /mi line. In addition, this line is 
very sensitive to the depth at which the dust becomes optically 
thick and a quantitative comparison to the far-infrared lines is 
thus difficult. 

2.1. Parameters for the HD 100546 model 

To model HD 1005 46, we chose t o use the dust density distri- 
bution obtained bv lMulders et al.l (1201 ll) . They fit the SED us- 
ing a dust radiative transfer calculation with a vertical structure 
that obeys hydrostatic equilibrium assuming that r gas = Td ust . 
Dropping this assumption results in a density structure w ith a 
puffed-up inner rim with about 10 AU dKamp et alJl2010h . but 
this region does not contribute significantly to the emission we 
probe here. For consistency, we also use the dust temperature of 



Mulders e t al.1 (1201 ll) . The surface density power-law is taken to 
be 2 oc r , which fits the VISIR images. Anothe r dust density 
struct ure for the HD 100546 disk was proposed by Benistv et al. 
(120101) . based on a analytical disk structure and fitted to the IR 
interferometric data. Their ma in conclusions abo ut the inner disk 
were incorporated into the iMulders et al.l d201 lb model that we 
use. It is well-known that SED fitting yields degenerate solu- 
tions and the obtained dust structure does not necessarily reflect 
the gas structure, e.g. owing to dust settling effects. In addition, 
the st r ucture of the disk i s likely non-axisymmetric (Pani c et al. 
| 2010t lOuanz et all 1201 lb and with substructures (Ardilaetal. 
2007). However, this approximate approach is a reasonable one 
for studying the main characteristics of the gas emission, given 
the large uncertainties in the gas line modeling (Section l53i l. 

We adopted a Hipparco s distance of 103 pc to HD 100546 
(Ivan den Ancker et al.ll 199 8). The inclinatio n and position angl e 
were taken to be 42° and 145°, respectively dPantin et al.l l2000). 
We used a Ke plerian velocity fiel d around a 2.5 M Q star, which 
was found by Panic et al] d2010h to closely reproduce the ob- 
served spectra. The outer radius is uncertain but estimated to 
be be tween 200 AU and 400 AU dPantin et al.ll2000t iPanic et al.l 
120101) . 

Inside the large inner hole of the HD 100546 disk, a 
smaller inner disk e xtending from 0.25 AU to 4 AU was found 
dBenistv et al.l l2010). We run models including the inner disk, 
but found that more than 99 % of the flux of the lines discussed 
in this work emerges from the outer disk. We thus considered 
only the outer disk starting at 13 AU for this work. The dust 
opacity of the inner disk was however included to calculate the 
local UV field. 

As input stellar radiation field (FigureQ]), we analyzed dered- 
dened FUSE observations for wa velengths 911-1 200 A, which 
had been previ ously analyzed by iDeleuil et al.l d2004l) and IUE 
observations dValenti et a l. 2003) for longer wavelengths (1200- 
3200 A). T he spectra were de reddened using the extin ction of 
Ay = 0.36 (lArdila et al.ll2007l) and the extinction law of lDraind 
d2003l) for ^ v = 3.1. The FUSE/IUE spectra were extended t o 
longer wavelengths using the B9V template of lPicklesI d!998l) . 
The X-ray luminosity h ad been found to be log(Lx) = 28.9 erg 
s (IStelzer et al.ll2006l) . Since the plasma temperature had not 
been well-constrained by observations, we assumed a rather high 
temperature of kTx = 6 keV. This hard spectra allows X-rays to 
penetrate deep into the atmosphere. Nevertheless, the lines stud- 
ied in this work change by less than ~ 5 % by switching on 
X-rays. Thus, the assumption on the shape (plasma temperature) 
of the X-ray spectra does not affect our results (Section 13.1b . 
In addition to the stellar r adiation field, we also account for the 
interstellar radiation field (lDrainelll978l) . Another source of ion- 
ization are cosmic rays, for which a rate of £ c . r . = 5 x 10~ 17 s _1 
is adopted. 

The dust opacities control the penetration of the photons and 
thus both the photodissociation and ionization of molecules as 
well as t he dust temperature. The dust grains used in the SED 
fitting by Muld ers et al.1 (1201 ll) are silicates with a 5 % (in mass) 
composition of carbonaceous material. For gas/dust = 100 and 
a present-day solar p hotosphere carbon abundance of 2.7 x 10 4 
relative to hydrogen dAsplund et al.l2009l) . this amounts to about 
15 % carbon bound in this population of grains. The grains 
have a size of 0.1 - 1.5 fim with a distribution proportional to 
a~ 3 - 5 , e mploying the Ma this, Rumpl and Nordsieck size distri- 
bution dMathis et al.lll977l) . The optical constants are for irreg- 
ularly shaped distributions of hollow spheres (DHS, iMin et all 
2005). The DHS dust model has been adopted to be consistent 
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Fig. 1. UV spectra of HD 100546 at the source distance of 103 
pc. We indicate the far ultraviolet (FUV) wavelength range im- 
portant for gas heating and the wavelength range of photodisso- 
ciation (p.d.) and photoionization (p.i.) cross-sections for differ- 
ent molecules. 

with the SED fitting of lMulders et alj d2011l) . This dust popu- 
lation has a mass of IP" 4 Af n , w hich agrees with observations 
of small grains bv lDominik et alj d2003l) . However, much of the 
mass can be in larger grains, and the gas/dust ratios given rela- 
tive to this population of grains are thus upper limits. This pop- 
ulation of grains does not contain the UV/visual attenuation by 
very small grains (VSGs) and polycyclic aromatic hydroca rbons 
(PAHs ). We thus also use R\r = 3.1 and 5.5 opacities after Draine 
(120031) at wavelengths < 1 fim (Figure to include VSG/PAH 
absorption. Using these opacities at short wavelengths only af- 
fects (i.e. increases) the dust temperature in the uppermost layer, 
above the region from which the line emission discussed in this 
work emerges, as we checked. For simplicity, w e thus use the 
same dust temperature for all models taken from lMulders et af] 
J201ll) . 




0.1 0.2 0.5 1 2 5 10 20 50 100 
Wavelength (^tm) 



Fig. 2. Adopted dust opacities. The solid lines give absorption 
mass extinction coefficients, while dashed lines indicated scat- 
tering mass extinction coefficients. 

The PAH abundance affects both the heating by means of the 
photoelectric effect, as well as the ionization balance, by means 



of recombination reactions. We adopt a 5 % PAH-to-dust ratio 
within 100 AU and 1% outside 100 AU, following a simulta- 
neous fitting of Spitzer an d VISIR N-band ob servations using 
the N c = 100 opacities of iDraine & "□ d2007l) (Gijs Mulders, 
pers. comm.). The PAH-to-dust mass ratio of 5 % corresponds 
approximately to the average Galactic PAH-to-dust mass ratio 
of 4.6 %, which in turn represents about 5 ppm of C per H 
nucleon bound in PAHs dDraine & Lil 120071) and 15 % of the 
present-day solar photosphere carbon abundance (for gas/dust 
= 100). The PA Hs are only excited by UV in the uppermost 
layer of the disk ( Visse r et al.ll2007l) and the observed features 
thus trace this region. In addition to the large uncertainties in the 
molecular properties of the PAHs, the PAH abundance is only 
poorly determined in the region from which the lines discussed 
in this study emerge. We thus choose to scale the PAH abun- 
dance with the gas/dust ratio, such that gas with gas/dust = 20 
contain s five times smaller amount of PAHs than for gas/dust 
= 100 dJonkheid et alJ l2007). We examine the influence of the 
PAH abundance separately in Section 14.31 We assume PAHs 
to be photodestroyed in regions with integrated far ultraviolet 
(FUV) fluxes > 10 6 ISR F, following t he resu lts of a more com- 
plete PAH chemistry by IVisser etal] d2007l) . This affects how- 
ever only the upper layer of the disk and we find that the lines 
discussed in this work are unaffected by this assumption. We 
note that the emission of all lines studied here emerges from zjr 
lower than the region of PAH destruction (Section [3~Tll3.4t . 



2.2. Observations 



Spectrally resolved observations of the low-/ CO 7 = 7-6, 
7 = 6-5, and 7 = 3-2 lin e s wer e obtained with the 12m 
APEX telescope by [Panic et al.l (1201 Oh and displayed a double- 
peaked line shape, that is consistent with a rotating disk. We note 
that the uncertainty in the 7 = 7-6 line is considerable and 
the 7 = 6-5 line is detected at a much higher signal-to-noise 
ratio. They have also searched for the 370 fim (809 GHz) [C I] 
3 P2 - 3 7 > i line, but were unable to detect it. We note that they 
have been using position switching with an off position 16' from 
the source. Thus, it is unlikely that the emission at the source 
had been canceled out by extended emission. The 609 /jm (492 
GHz) 3 7'i - 3 Po line of neutral carbon was not observed. 

The [CII] 158 ^m 2 P 3/2 - 2 Pi/ 2 , the [OI] 63 /mi 3 P X - 3 P 2 , 
and 145 yum 3 Po - 3 P; lines and mid-7/hig h-7 l i nes o f CO with 
7 up = 14 to 30 were observed by ISturm etall d2010t) with the 
PACS instrument onboard Herschel as part of the "Dust, ice, 
and gas in time" (DIGIT) key program (PI: N. Evans). The CO 
7 = 31 - 30 line was blended with the OH 3/2,7/2+ - 5/2- 
line and CO 7 = 23 - 22 was potentially blended with H2O 
4i4 - 3o3- These two lines are thus taken as upper limits. We 
note that the apparent jump in CO line flux at 7 > 25 may have 
been introduced by uncertainties in the PACS flux calibration 
when the first PAC S results were published, including those of 
ISturm etakl (|2010). The larger uncertainty in these lines does 
however not affect any of our conclusions. The PACS data are 
spectrally unresolved. Except for [CII], all lines are consistent 
with their emergence from within the central 9'.' 4 pixel. To cor- 
rect for extended emission in [CII], we subtracted the pixels 
neighboring the central pixel. We also accounted for emission 
leaking from the center pixel owing to the extended point spread 
function. The origin of the [C II] line is discussed in Section l331 
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Table 1. Line fluxes measured towards HD 100546. 



Species 


Line 


Wavelength" 
[pm] 


Flux/Error'' 

1 n— 1 7 rii t — 9 n 

10 LW m -J 


APEX (Panic et al. 2010) 






CO 


7 = 3-2 


866.96 


0.0149 ± 0.003 


CO 


7 = 6-5 


433.56 


0.132 ± 0.007 


CO 


7 = 7-6 


371.65 


0.11 ± 0.02 


[CI] 


3 P 2 - 3 P I 


370.42 


<0.009 


Herschel PACS (Sturm et al. 2010) 




CO 


7 = 14- 13 


186.00 


7.4 ± 0.9 


CO 


7= 15-14 


173.63 


11.5 ± 0.8 


CO 


7= 16- 15 


162.81 


8.3 ± 0.9 


CO 


7 = 17- 16 


153.27 


10.5 ± 0.8 


CO 


7 = 18- 17 


144.78 


9.9 ± 1.0 


CO 


7 = 19- 18 


137.20 


8.9 ± 0.7 


CO 


7 = 20 - 19 


130.37 


7.1 ± 0.7 


CO 


7 = 21 - 20 


124.19 


5.9 ± 0.6 


CO 


7 = 22 - 21 


118.58 


5.9 ± 0.9 


CO 


7 = 23 - 22 


113.46 


<4.8 ± 1.2 


CO 


7 = 24 - 23 


108.76 


10.6 ± 1.3 


CO 


7 = 26 - 25 


100.46 


12.7 ± 2.4 


CO 


7 = 27 - 26 


96.77 


1 n O a. 1 A 


CO 


7 = 28 - 27 


93.35 


10.0 ± 1.0 


CO 


7 = 29 - 28 


90.16 


11.8 ± 1.6 


CO 


7 = 30 - 29 


87.19 


8.5 ± 1.0 


CO 


7 = 31-30 


84.41 


<4.9 ± 1.5 


[cm 


2 Pl/2 — 2 Pl/2 


157.74 


13.5 ± 1.5 C 


[OI] 


'Pi-'Pl 


63.18 


554 ±5 


[OI] 


3 Po- 3 P l 


145.53 


35.7 ± 1.3 



Notes. <a) Rest frequency ( ' Formal errors, not including systematic er- 
rors <r) Extended emission subtracted, see text. 



3. Results: A representative model 

We now present the details of a representative model of HD 
100546. By representative, we mean that it reproduces suffi- 
ciently well the observations enabling us to study the main phys- 
ical and chemical effects. It is however not meant to be a best fit 
model. The dependence on the parameters assumed here is dis- 
cussed in the next section. The outer radius of the representative 
model is assumed to be 400 AU. We define the depletion fac- 
tor 6c as the carbon abundance in the gas phase relative to the 
present-day s olar photosphere abu ndance of 2.7 x 10 relative 
to hydrogen dAsplund et al.ll2009l) . In this section, we assume a 
carbon gas phase abundance of 1.3 x 10~ 4 relative to hydrogen 
(6c ~ 0.5), which is characteristic of the diffuse ISM. Adopting 
gas/dust = 20, we get a gas mass of ~ 2 x 10~ 3 M©. We use the 
Rv = 5.5 dust extinction law for this model. 

3.1. Physical structure 

The physical conditions of the representative model are shown in 
Figure[3] Only one quadrant of the disk is given and only regions 
with gas density > 10 5 cirT 3 are shown. We present the gas den- 
sity in Figure [3^/b both in linear (r,z) as well as in (\og(r),z/r)- 
space. The latter has the advantage that the inner part and upper 
disk can be seen clearly. Direct rays from the star to the disk cor- 
respond to vertical lines. The density reaches up to a few times 
10 10 cm~ 3 in this model in the inner part of the disk midplane. 

The X-ray ionization rate is given in Figure [3p. Despite the 
low X-ray luminosity of this source of ~ 8 x 10 28 erg s , the 
ionization rate due to X-rays in the inner disk is still about three 
orders of magnitude higher than the standard cosmic ray ioniza- 
tion rate of 5 X 10~ 17 s _1 . For the heating and cooling budget of 



the disk, X-rays can however be neglected, since the luminos- 
ity in the FUV band (6-13.6 eV) is at ~ 4 x 10 34 erg s" 1 much 
higher than the X-ray luminosity. Nevertheless, the considerable 
ionization rate may affect the chemistry and for example lead 
to a destruction of water (e.g. IStauber et alj|2006l) . The simple 
species considered in this work are however little affected by X- 
rays and we find that the CO, C, C + , and O fluxes change by less 
than ~ 5 % when X-rays are switched off. 

The (i ntegrated) FU V flux in units of the interstellar radia- 
tion field (Drainel ll978t 1 ISRF ~ 2.7 x 10 3 erg s 1 cm 2 ) is 
shown in Figure[3]i Very high FUV fluxes of up to 10 8 ISRF are 
reached in the surface of the inner edge of the disk. In the upper 
atmosphere of the outer disk at r = 400 AU and z/r — 0.3 the 
FUV flux still reaches 10 5 ISRF at a density of a few times 10 5 
cirT 3 . We note that these are rough l y the conditions of the Orion 
bar PDR dHogerheiide et al.l [19951 ijansen et alJll995h . Similar 
conditions are also predicted for the outflow walls, the "dense 
PDR" along the outfl ow of a high-mass star forming region 
dBruderer et alJl2009al) . 

Figure^ gives the FUV extinction in units of Ay, obtained 
by defining tfuv = ~ log(/attMinatt)> with the attenuated FUV 
flux 7 att and th e FUV flux corrected only for geometrical dilu- 
tion / unatt (e.g. lBruderer et all2009aL Eq. 2). Conversion to Ay is 
done by scaling with the dust opacities at 2070 A and 555 A and 
accounting for Ay ~ 1.086rv. The UV extinction is not used 
during the calculation, because the photoionization and disso- 
ciation rates are calculated from the wavelength-dependent UV 
spectra. The FUV extinction however shows that owing to scat- 
tering, FUV photons can penetra te much further into the d isk, 
particularly in the outer parts (e.g. lvan Zadelhoffet al1 l2003h . 

Figure [3f shows the adopted PAH abundance. The drop in 
the inner, upper disk is due to the photodissociation of PAHs. 
The drop at 100 AU is observationally constrained (Section lXTT i. 

The dust and gas temperature are presented in Figure [3g/h. 
Figure[3j shows the ratio of the gas to dust temperature. The dust 
temperature does not exceed ~ 250 K, but stays above ~ 30 K 
across the entire model. Thus, CO does not freeze out onto the 
dust grains. The gas temperature exceeds the dust temperature 
by factors of up to 50 in the inner upper disk and reaches tem- 
peratures of up to ~ 6000 K. In the outer upper disk, the gas 
temperature is still r gas ~ 1000 K and thus about a factor of ten 
higher than 7d us t- Deeper within the disk, below z/r — 0.15, col- 
lisions between gas and dust grains result in a coupling of the 
gas and dust temperature (T gas = 7d us t)- 

3.2. Chemical structure 

Fractional abundances of C + , C, CO, O, H2, and electrons rel- 
ative to the total hydrogen density («h = "(H) + 2n(H2)) are 
shown in Figure [4] As found by previous chemical models of 
protoplanetary disks the chemical structure resembles a classical 
PDR in the vertical direction (e.g. iJonkheid et al.|[2007h . Below 
a top layer of atomic hydrogen, H2 forms, as self-shielding be- 
comes sufficient. In the outer disk, the inward column is also 
large enough to prevent H2 from photodissociation and hydro- 
gen is in molecular form. 

Carbon in the upper disk (z/r ~ 0.2) is mainly in the form of 
C + , followed by neutral carbon and CO deeper within the disk. 
An exception is the "warm finger" of CO in the outer disk (r > 
100 AU, z/r > 0.25). In this region, CO forms through the re- 
action of C + with H2 to CH + , followed by a reaction with O to 
form CO + which t hen reacts with H2 t o form HCO + and sub- 
sequently CO (see lJonkheid et ai1l2007l) . The initiating reaction 
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r (AU) r (AU) t (AU) 

Fig. 3. Physical conditions in the representative model. The contour lines are shown for values indicated in the label by small arrows. 
The thin dashed line indicates z/r = 0.15,0.2, and 0.25. The panels show: a.) Gas density b.) Gas density in linear scale, c.) Total 
ionization rate, d.) Local FUV flux e.) FUV extinction f.) PAH abundance in units of the ISM abundance g.) Dust temperature h.) 
Gas temperature i.) Ratio of the gas to dust temperature. 



is strongly endothermic with an energy barrier corresponding to 
4640 K, thus requires a high temperature to proceed. The OH 
formation at high temperature, by means of the reaction of O 
and H2, leads to additional CO + for this chain of reactions. The 
OH reacts with C + to form CO + . At the inner edge of the disk at 
13 AU below z/r ~ 0.1, CO is shielded from photodissociation 
owing to shadowing by the inner disk. The higher abundance of 
CO at r = 13 AU and z/r — 0.2 is due to the adopted H2 forma- 
tio n rate: since the stic king coefficient decreases above ~ 3000 
K (iCuppen et alJl20 lOf) . the lower gas temperature at higher z/r 
leads to an increase in H2 and subsequently in CO. 

Atomic oxygen has a considerably high abundance through- 
out the disk, except for the midplane. This is because the elemen- 
tal abundance of oxygen is higher than that of carbon and not all 
oxygen can be locked into CO, which is almost unaffected by 
photodissociation, by means of self-shielding. In addition, oxy- 
gen cannot be ionized directly by FUV photons since its ioniza- 
tion energy is slightly higher than 13.6 eV. In the inner disk with 



13 AU < r < 20 AU and z/r < 0.1, gas and dust are at a tempera- 
ture > 100 K and the gas is not exposed to strong FUV radiation. 
In this part of the disk, oxygen is locked into water. Further out, 
water freezes out and locks most of the oxygen into water ice. In 
this region, the carbon is in the form of gaseous CH4. 



The electron fraction roughly follows the abundance of C + . 
Charge exchange with PAHs is another important reaction type 
for the ionization balance. The PAHs are positively charged in 
the upper atmosphere, neutral at intermediate heights, and nega- 
tively charged at greater depths within the disk due to charging 
by electron attachment. Since the electron abundance is very low 
in the midplane, the PAHs can again become neutral. Charge ex- 
change C + + PAH -> C + PAH + and recombination C + + PAH" 
— > C + PAH are the two reactions that affect the simple carbon 
chemistry most, leading a higher abundance of neutral carbon. 
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r (AU) r (AU) r (AU) 

Fig. 4. Fractional abundances of C + , C, CO, O, H2, and electrons relative to the total hydrogen density («h = "(H) + 2n(H2)). The 
thin dashed line indicate z/r = 0.15, 0.2, and 0.25. Abrupt changes at R — 100 AU are due to a varying PAH abundance, see text. 



3.3. Warm atmosphere versus dust temperature 

Line fluxes predicted by the representative model are presented 
in Figure [5] The figure gives the CO ladder on the left and the 
[CI], [CII], and [OI] atomic fine-structure lines on the right. 
Atomic fine-structure lines of C and CO lines up to J — 7 - 6 are 
convolved to the APEX beam and all other lines to the Herschel 
beam. In addition to the model with the calculated gas tempera- 
ture, a model that assumes that the gas temperature is equal to the 
dust temperature is shown. A factor of two around the observed 
fluxes is shown by a gray shaded region. We consider this region 
as a good agreement to observations given all the uncertainties 
in data and modeling (see Section lBTBI l. 

The model with the calculated gas temperature fits the data 
well, with the exception of the high-/ CO lines with J > 25 
(see Sections 14.51 and l33T l. However, the cool atmosphere with 
r gas set to Jdust fails to reproduce the high-/ CO ladder. The 
/ = 16-15 transition with E up ~ 750 K (Tabled is already un- 
derproduced by about an order of magnitude, and the higher-/ 
lines are underpredicted even more. For example, / = 25 - 24 
(£ up ~ 1800 K) is underpredicted by more than two orders of 
magnitude. The atomic fine-structure lines have upper level en- 
ergies that are much lower than the high-/ lines of CO. Both the 
[CI] and [CII] lines have E up < 100 K and the line fluxes cal- 
culated with a cool and warm atmosphere are very similar. The 
[C II] line is somewhat underproduced and [C I] overproduced. 
However, [O I] has an upper level energy of 327 K (145 fim) and 
227 K (63 /mi). Consequently, the [OI] lines for the cool atmo- 
sphere differ by an order of magnitude from those obtained for 
the warm atmosphere. We conclude that the prominent emission 
in high-/ CO and [O I] atomic fine-structure lines is strong evi- 
dence of a higher temperature of the gas in the upper at mosphere 
of the disk. This agrees to the recent detection of CH + dThi et al.l 



120111) and OH dSturm et al.ll2010l) towards this disk, which are 
both formed in endothermic reactions with H2. 

3.4. Origin of the line emission 

Where does the observed line emission originate? The angular 
resolution of Herschel, corresponding to about 1000 AU for 63 
fim at the distance of HD 100546, does not resolve the disk. The 
model however, may help us to locate the region of emission. 
Figure [6] gives the relative contribution to the observed fluxes. 
We show the integrated contribution function (ICF), which in- 
corporates all the effects involved in the line formation: abun- 
dance of the molecule, excitation, and opacity. In addition, we 
provide the two radii delineating the region from which 75 % of 

the emission emerges, either inside or outside. 

The ICF, disc u ssed in iTafalla et alj d2006l) and 
Pavl vuchenkov et al] J2008), shows the relative contribu- 
tion to the line emission along a ray. We obtain from the formal 
solution of the radiative transfer equation (Eq. IA.14I ) 

ICF= f (S v -S dusl )e^ (l-e-^)dv, (1) 

«^line 

where the source function is S v , the dust source function is S dust> 
the opacity to the observer is r v and the opacity within one cell is 
Ar v . The integration is performed over frequency. For simplicity, 
we take the disk to be face on and allow the rays to propagate 
parallel to the z axis (as e.g. in iPanic et al.l 120 10b . This avoids 
complications with the adopted Keplerian velocity profile but 
conserves the main physical conclusions. The contribution to the 
observed line flux from an annulus located at a radial distance r 
is given by 2nl(r)rdr, with the flux I{r) at radius r. Thus, the rel- 
ative contribution to the observed flux is oc I(r)dr. In Figure [6] 
we thus show the contribution function times the radius (rxICF), 
scaled to its peak value. 
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Fig. 5. Integrated line fluxes obtained 
from the representative model. The left 
part of the figure shows the CO ladder 
and the right part atomic fine-structure 
lines. Observations are indicated by 
black crosses. The red line and squares 
show model fluxes for a model with cal- 
culated gas temperature, the blue line 
and triangles with r g!ls set to Tdust- The 
gray shaded region indicates a factor of 
two to the observed fluxes. 



The low-/ emission of CO 7 = 3-2 originates at radifl 
between r ~ 70 - 220 AU at z/r ~ 0.15 - 0.25. The density 
above this region is too low to contribute significantly to the 
emission, since the regions below are shielded by the line opac- 
ity. Thus, the far side of the disk does not contribute much to 
the emission, except for some radiation in line wings. Owing 
to the increasing density and temperature towards the inner re- 
gion at a given z/r, even spatially smaller regions at radii < 100 
AU contribute to the observable emission. Mid-7 lines, e.g. CO 
J = 16-15 with E up ~ 750 K, are not excited in the outer disk 
and only molecules within r ~ 35-80 AU contribute to the emis- 
sion. These lines reach line opacities of around unity and the far 
side of the disk does not contribute to the emission in the inner- 
most regions (radii < 30 AU). The highest-/ emission, e.g. CO 
J = 30 - 29 with £ up ~ 2800 K, detected towards HD 100546 
emerges from the inner r ~ 20 - 50 AU in a thin layer at the 
surface of the disk. There, CO is formed but the gas temperature 
is still sufficiently high to excite the lines. The opacity of these 
lines is optically thin and the dust extinction is not large. Thus, 
emission from both the near and far sides reaches the observer. 

The atomic fine-structure lines of [CI] and [CII] emerge 
mainly from the outer disk at radii ~ 150 - 300 AU. Owing to 
the higher abundance of those species in the upper atmosphere, 
combined with the low critical density (~ 10 3 cm 4 ) and low up- 
per level energy of the [C I] and [C II] lines, these lines emerge 
from regions higher in the disk than the low-/ lines of CO. The 
oxygen [O I] lines have considerably higher critical densities of 
~ 10 5 - 10 6 cm -3 and also higher upper level energies. They are 
thus not excited in the outer disk and the line emission comes 
from regions within ~ 50-210 AU. Since this line is optically 
thick, most of the emission is from the near side of the disk. 

3.5. CO line profile 

The predicted shapes of the CO / = 3-2,/ = 10-9,/ = 16-15 
and / = 30 - 29 lines are shown in Figure [7] For CO / = 3 - 2, 
the profile is given for the APEX beam, while the higher-/ lines 
are calculated for the Hersche l beam. The HIFI hetero dyne spec- 
trometer onboard Herschel (ide Graauw et al.1 [201 Oh can spec- 
trally resolve lines down to < 0.1 km s _1 and covers the fre- 
quency range of CO / = 5 - 4 up to CO / = 16 - 15. For 



1 Measured by the 75 % contribution radii given in Figure[6]and Table 

® 



Table 2. Line properties and origins 



Line 




A 


P 

£ -*up 


Aui 


Wcrit" 


Origin 






1pm] 


[K] 


[s- 1 ] 


[cm- 3 ] 


[AU] 


CO / = 3 - 


2 


867 


33 


3(-6) 


2(4) 


70-220 


CO / = 16 


- 15 


163 


752 


4(-4) 


1(6) 


35-80 


CO / = 30 


-29 


87 


2565 


2(-3) 


4(6) 


20-50 


[CI] Vi- 


3 P 2 


370 


62 


3(-7) 


1(3) 


150-300 


tC II] 2 Py 2 


- 2 P 1/2 


158 


91 


2(-6) 


5(3) 


150-300 


[OI] 3 P, - 


3 P, 


63 


228 


9(-5) 


5(5) 


50-180 


[OI] 3 P - 


'Pi 


145 


327 


2(-5) 


5(4) 


80-210 



Notes. a(b) means a x 10*. (a) Critical density (A„//2/C„/) for colli- 
sion with H 2 at the temperature of E up . Atomic/molecular data from 
ISch6ieretalJd2005h . 



illustration purposes, we however also show CO / = 30 - 29, 
which might eventually be accessible to resolved observations 
with the GREAT instrument onboard SOFIA. 



— J = 3-2 
J = 10 - 9 




-15 -10 -5 5 10 15 

Velocity (km s _1 ) 

Fig. 7. Predicted line shapes for the CO / = 10 - 9, / = 16 - 15, 
and J - 30 - 29 in the beam of Herschel. The CO / — 3 - 2 line 
predicted for APE X is shown together with the observations of 
lPanicetal.1 (120101) . 
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Fig. 6. Contribution function to the line formation. The disk is viewed from the top and the contribution function is normalized to 
the peak. Total (line+dust) opacities of Tij necent er = 1 at the line center are indicated by a black line, and dust opacities of Tdust = 0.1 
by gray lines. The red arrows indicate the radii from which 75 % of the emission emerges, either inside or outside. The thin dashed 
line indicates \z/r\ = 0.15,0.2, and 0.25. 



The modeled CO 7 = 3-2 l ine agrees well in terms of 
width with the spectra observed bv lPanic et al.l (l2010h . A slight 
asymmetry in the observed spectra with a st ronger blue-shifted 
peak is however not reproduced by the model. [Panic et al.l (2010) 
attributed the asymmetry to an asymmetric te mperature struc - 
ture that might be due to a warped inner disk (lOuillenl 12006). 
Modeling these asymmetries is beyond the scope of this study 
and would require tighter observational constraints on the den- 
sity structure of the disk. 

The CO 7 = 16 - 15 and 7 = 30 - 29 lines are predicted 
to be considerably broader than those of CO 7 = 3-2 and 
j = 10-9. While CO 7 = 3 - 2 and 7 = 10 - 9 both have 
a width of ~ 5 km s _1 , CO 7 = 16 - 15 and 7 = 30 - 29 
have widths of ~ 9 km s _1 and ~ 12 km s _1 . This can be un- 
derstood based on the origin of those lines from regions closer 
to the star. For the assumed Keplerian velocity field around a 
2.5 M star, the projected velocity along the semi-major axis is 



»proj 



sin(/) VM,G/r 



For a given velocity 



in the observed spectra, r corresponds to the maximum ra- 
dius from which the line photons can emerge. For the half width 
of the CO lines discussed here (2.5, 4, and 6 km s -1 ) we get 
r = 160, 50, and 27 AU, roughly corresponding to the radii of 



the main emission found in Section 13.41 We note that the CO 
7 = 30 - 29 line is broader than expected from the projected 
Kepler velocity at the inner rim. This is due to thermal broad- 
ening of the line, corresponding to more than 1 km s _1 for the 
gas temperatures of a few 1000 K reached in the inner, upper 
atmosphere. 

We conclude that the width of the CO lines is an interesting 
indicator of the radial origin of the emission and thus a probe/test 
of the temperature structure. With current facilities, lines up to 
CO 7 = 16-15 can be spectrally resolved, allowing us to study 
regions within radii of a few tens of AU as suggested by our 
results. 



4. Results: Dependence on parameters 

This section presents the results of a parameter study. We vary 
different input parameters of the model, such as the size of 
the disk or the amount of gas in the disk in order to under- 
stand the dependence of the line fluxes on those parameters. 
Understanding these dependences will help us to interpret obser- 
vational trends and also to quantify the robustness of the results, 
given that large uncertainties enter the modeling (Section lBTBT ). 
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4.1. Varying size, gas/dust ratio, and carbon abundance 

Key parameters for the gas modeling are the total amount of gas 
and the size of the disk. In this section, we vary the outer ra- 
dius of the gas disk to be either 200 or 400 AU and assume that 
gas/dust = 20 or 100. A disk outer radius of r out = 400 AU 
and gas/dust = 20 (= 100) yields a gas mass of 8 x 10~ 4 M 
(4 x 10~ 3 M Q ). For an outer radius of r out = 200 AU, the gas 
mass is 4 x 10~ 4 M Q (2 X 10~ 3 M ), thus a factor of two lower 
than for r out = 400 AU. 

The o uter radius of the d isk is difficult to determine, as dis- 
cussed in Mulders etal.1 (1201 lh . For H P 100546, scattere d light 
has been detected out to ~ 1000 AU dArdila et al.ll2007l) . but it 
remains unclear whether this emission originates from the disk 
or a diffuse remnant envelope. Thus, the dis k might be s ignifi- 
cantly smaller. Modeling the low-/ CO lines. iPanic et al. 
constrained the size o f the gas disk to be ~ 400 AU and the gas 
mass to be 10 -3 M . iMulders et al.l (1201 lh employed an expo- 
nential cutoff in sur face density outside 350 AU. In other disks, 
Hug hes et al.l d2008) found that CO J = 3 -2 can be more clearly 
explained when such a exponenti al cutoff is used. Here, we use 
a sharp cutoff at r out = 400 AU dPanic et al.l | 2010l). but we als o 
consider a disk with a radius of only 200 AU dPantin et al.l2 000). 

The gas/dust ratio is assumed to be either an ISM value of 
100 or a lower value, as previously suggested by Cha pillon et al.l 
(1201 Oh to explain the non-detection of [C I] 370 /urn and 609 /im 
lines towards the Herbig Ae star CQ Tau. We adopt gas/dust = 
20, which is similar to their choice. As mentioned before, these 
gas/dust ratios only refer to the population of small grains and 
are thus upper limits. 

The carbon budget of the gas depends on the elemental de- 
pletion on dust grains. The warm ISM shows signs of little or 
no depletion of oxygen and we thus assume that all oxygen not 
locked up in silicates is in the gas phase. We note that oxygen 
may freeze-out as molecules onto the dust grains, e.g. as water 
ice, but this happens only in the cold shielded regions of the disk 
(see Section [3~TT i. Carbon, however, may be depleted by a consid- 
erable fraction by locking it up in refractory material. We chose 
a carbon abundance of 6c ~ 0.1, corresponding to an abundance 
of 2.4 x 10~ 5 relative to hydrogeifl In this way, we can study the 
degeneracy between having a higher gas/dust ratio of 100 com- 
pared to 20 and a lower carbon gas phase abundance of <5c =0.1 
compared to 6c ~ 0.5. For gas/dust = 100, we also run models 
with 6 C = 0.05. 

In Figure [8] we show the integrated line fluxes for the mod- 
els with different radii, gas/dust ratios, and carbon abundances in 
the gas, as in Figure [5] We first discuss the models with gas/dust 
= 20 shown in the upper panel of the figure. The low-/ and high- 
/ lines of CO have different behaviors. As expected, the high-/ 
lines do not depend on the outer radius of the disk since the radi- 
ation emerges from the inner disk. The optically thin high-/ lines 
however scale with the amount of carbon in the gas. We note that 
this is already true for mid-/ lines that are marginally optically 
thick. The optically thick low-/ lines depend less on the amount 
of carbon, with a change between the model of the same radius 
of about a factor of two. We note that these optically thick lines 
may still depend on the abundance, since e.g. the line wings may 
remain thin. The low-/ lines also depend on the outer radius of 
the disk. Changing the radius from 200 AU to 400 AU increases 
the area of the disk by about a factor of four. The low-/ lines 
however only change by about a factor of two. This is because a 



2 Sq is defined as the fraction of carbon in the gas-phase (See Section 

13 



considerable part of the emission originates from within 220 AU 
(SectionED. 

As expected, the fine-structure lines of [OI] for gas/dust 
= 20 do not change with the amount of carbon in the gas. They 
however slightly scale with the outer radius by about the same 
factor as the low-/ lines of CO. This can be understood by these 
lines emerging from within about the same radial region, as dis- 
cussed in Section [3~4l 

The [C I] and [C II] lines depend on both the outer radius and 
the amount of carbon in the gas phase by similar amounts. This 
is the result of both lines being optically thin and emerging from 
approximately the same region. The line fluxes indeed scale with 
the amount of carbon in the atmosphere, while they change by 
about a factor of four with radius. 

The models with gas/dust = 100, shown in the lower panel 
of Figure[8] exhibit very similar dependences on the outer radius 
and amount of carbon in the gas phase. How do the lines change 
with the larger amount of gas in the atmosphere? The optically 
thin high-/ lines of CO approximately scale with the amount of 
gas. However, since the thermal balance is also affected by the 
higher density, leading to slightly higher temperatures, the lines 
increase by more than a factor of five. As expected, the low-/ 
lines scale by a smaller factor, owing to optical depth effects. 
The [OI] fluxes approximately scale with gas mass. The carbon 
fine-structure lines scale by a smaller factor, as the higher density 
moves the C + /CO transition up in the atmosphere resulting in a 
smaller increase in the column density. 

The line ratio [C I]/[C II] is about a factor of two lower for the 
higher gas/dust ratio. This is due to different branching ratios in 
the C + destruction. The radiative association of C + with H2 leads 
to CHJ , while grain surface reactions of C + with hydrogenated 
PAHs (PAH:H) leads to CH + . Following a chain of reactions 
with H2 and electron recombinations, CH + and CHJ lead to CH 
and CH2. For higher densities, CH and C H? are turned into CO 
rather than being photodissociated (see iTielens & Hollen bach 
1985). This only slightly shifts the CO transition, thus affects 
the CO fluxes only slightly, but does affect the [C I] lines. 

There is some degeneracy between models having more gas 
(gas/dust = 100) but less carbon (6c = 0.1) and models with less 
gas (gas/dust = 20) and more carbon [6c = 0.5), since they have 
the same amount of carbon in the atmosphere. These models 
with the same outer radius agree more closely than about a fac- 
tor of two in the low-/ to mid-/ lines of CO and [C II] (see also 
Figure|9] models with Ry - 5.5). The high-/ lines of CO exhibit 
larger changes. As discussed before, neutral carbon changes by 
about a factor of three. 

4.2. Varying the dust opacity 

The dust opacities at UV wavelengths are another key parame- 
ter of the gas, as they control the penetration of the high-energy 
photons into the disk atmosphere. This affects both the chemistry 
and thermal balance, for example by means of the photodissoci- 
ation/ionization of molecules or the photoelectric effect on the 
dust grains. The opacity at short wavelengths is dominated by 
VSGs and PAHs. 

Observationally, the size of the smallest grains is difficult to 
determine. For HD 100546, indications of g rain process i ng an d 
larger minimum grain sizes were found by iArdila etafl ((2007). 
Bouwma n et al.l d2003l) suggested that the size distributions vary 
with radii, employing a population of smaller grains in the inner 
40 AU and larger grains in the outer disk. Using such a grain 
distribution, Benist v et alj d2010t) fit the SED using silicate dust 
grains and a surface layer of 0.05- 1 //m grains at 13-50 AU and 
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b.) gas/dust = 100 
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Fig. 8. Integrated line fluxes for models 
with different gas/dust ratios, outer ra- 
dius of the disk, and carbon abundances 
in the gas phase, a.) With gas/dust = 20 
b.) With gas/dust = 100. 
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Fig. 9. Same figure as Figure [8] but ei- 
ther assuming R\ = 3.1 orRy = 5.5 dust 
opacities (Figure [2}. The outer radius of 
the models is 400 AU. 



an outer disk consist ing of 1 fim - 1 cm size grains. In contrast 
Muld ers et al.l (1201 lh manage to reproduce the SED by assuming 
a single grain distribution of 0. 1 - 1 .5 jjm size, but using a verti- 
cal structure in hydrostatic equilibrium rather than a fixed scale 
height. In addition to radially dependent dust properties, dust 
settling m ay also introduce a ver tical stratification of dust prop- 
erties (e.g. D'Alessio et al. 2006), having a similar effect on the 
gas as a larger scale height of the gas st ructure and thus a h igher 
gas/dust ratio in the upper atmosphere (iJonkheid et alJ 2006). 



In this section, we study the influence of the dust opacity at 
UV wavelengths on the line fluxes. A detai led study of settling 
and radially dependent dust properties (e.g. Birn stiel et al.l feOlO) 
is beyond the scope of this study and we only globally change 
the dust properties at UV/Visual wavelengths. We assume the 
dust opacity laws given in Section |2~T1 (Fi gure l2b . These are ei- 
ther ISM-like grains with Ry = 3.1, corresponding to a mini- 
mum grain size of ~ 50 A, grains with Ry — 5.5 corresponding 
to some coagulation and grain growth and considerably larger 
grains with a minimum size of 0. 1 /mi. The models with a min- 
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Fig. 10. Same figure as Figure [8] 
but assuming 0.1-1 .5 fim dust grains 
<Mulders etaLll201lh . a.) With gas/dust 
= 10 b.) With gas/dust = 100. 



imum size of 0.1 yum do not contain any VSGs/PAHs and the 
opacity at UV/visual wavelengths is thus considerably smaller. 
Photodestruction o f VSGs to PAHs in the upper atmosphere 
(Bern e et al . 2009) leads to some additional vertical change in 
the UV opacity, even tho ugh their mass extinction coefficients in 
the UV are similar (e.g. iPontoppidan et al.l l2007. Figure 2). We 
thus consider this section as an approach to studying the influ- 
ence of different UV penetration lengths on the chemistry, inde- 
pendent of whether the absorption is due to VSGs or PAHs or 
the absence of both. 

In Figure [9] we show the line fluxes obtained using the Ry = 
3.1 and Ry = 5.5 dust opacities. The outer radius of the models 
is 400 AU, and we adopt either gas/dust = 100 and 6c =0.1 or 
gas/dust = 20 and 6c = 0.5. All fine-structure lines change by 
less than 50 % by using the two different opacity laws. The low- 
J lines of CO change similarly little, but the high-/ lines deviate 
by up to about a factor of two. The changes are very similar to 
those for all other combinations of outer radii, gas/dust ratio, and 
carbon abundance. 

The results obtained using the much larger 0.1-1.5 yum grains 
are given in Figure [TO] For gas/dust = 100, the line fluxes are 
much larger than the models with Ry = 5.5 dust opacities, with 
increases of up to an order of magnitude. Thus, the observed 
line fluxes can no longer be reproduced, not even with the lower 
carbon abundance. To roughly fit the observations, we thus show 
models with gas/dust = 10 in a second panel. The dependences 
of the models with large grains on the outer radius and carbon 
abundance is similar those of to the Ry =5.5 models. 



We conclude that the dust opacity, particularly that of the 
small grains governing the penetration of the UV radiation into 
the disk atmosphere, has major implications for the line fluxes. 
What is the reason for this strong dependence? In Figure[TT| we 
show vertical cuts through the disk in density, temperature, FUV 
field, and fractional abundance of CO and C + for the models 
with different dust opacity laws. Comparing the results obtained 
with Ry = 3.1 and Ry =5.5 shows that while the gas temper- 
atures in the upper atmosphere are similar, the transition to the 
cold atmosphere, where gas and dust temperature are coupled, 
occurs at greater depth. Similarly, the C + /CO transition occurs 
at greater depths within the disk. Both results are caused by the 
FUV field that can penetrate deeper into the disk. While the tem- 
perature and abundance structures are qualitatively similar for 
the Ry = 3.1 and Ry = 5.5 models, they are considerably dif- 
ferent for the larger grains: The gas temperature is still coupled 
to the dust temperature in the inner midplane (50 AU). Further 
out (at 250 AU), the gas temperature however no longer couples 
to the dust temperature. This is consistent with the much larger 
UV field in the midplane, which also explains the shift of the 
C + /CO transition deeper into the disk. Thus, larger grains allow 
UV radiation to penetrate more deeply into the disk to regions 
of higher density and also to heat the gas to higher temperatures. 
The combination of both leads to the considerable increase in 
line flux. 
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Fig. 11. Vertical cuts through the disk at a radius of 50 AU (top panels) and 250 AU (lower panels). Green lines show values obtained 
using the 0.1-1.5 /mi grains, blue and red lines indicate models calculated with Ry = 5.5 and Ry = 3.1, respectively. a.)/d.) Gas 
density, gas an dust temperature b.)/e.) FUV field in units of the interstellar radiation field c.)/f.) Abundance of CO and C + relative 
to the hydrogen density. 



4.3. Varying PAH abundance 

Given the uncertain abundance of PAHs and the difficulties in 
simultaneously reproducing the [C II] fluxes and [C I] upper lim- 
its, we explore the effect of the PAH abundance on the line 
fluxes in this section. The PAH abundance affects the models 
in two ways. First, the ionization structure is changed by PAHs 
by means of recombination reactions such as C + + PAH — > C 
+ PAH + , with PAH being neutral and PAH + positively charged 
PAHs. Second, the gas heating depends on the PAHs by means 
of the photoelectric effect on small PAHs. 

The emission of PAH can be direc tly observed by its bands 
in the infrared (e.g. iGeers et alJ 12007). This emission however 
only traces regions where PAHs are excited by UV p hotons. 
An abundance and excitation study by Viss er et all (120071) shows 
that in disks around Herbig Ae stars, most of the emission from 
intermediate-size PAHs emerge from the uppermost layer, above 
the Ty = 1 surface. Thus, PAHs do not trace down to regions 
where the C + /C transition is affected by the PAH charge ex- 
change. The PAH abundance is thus rather uncertain. 

Figure[T2]presents line fluxes obtained with models that con- 
tain different PAH abundances in the outer (r > 100 AU) and 
inner (r < 100 AU) part of the disk. The outer radius of the 
disk is 400 AU, an ISM carbon abundance, and gas/dust = 20, 
corresponding to the representative model presented in Section 

m 



Changing the PAH abundance at radii > 100 AU from 20 % 
ISM to no PAHs, hardly changes the CO ladder. In addition, the 
[O I] fine-structure lines are hardly affected. The [C II] line how- 
ever slightly increases and [C I] decreases by a factor of about 
two, owing to the missing charge exchange converting C + to 
C. Changing the PAH abundance in the outer disk to the ISM 
abundance increases the mid-7 lines of CO, owing to the greater 
heating caused by the higher abundance of PAHs. We note that 
the high-/ lines of CO are less affected than the lower-/ lines, as 
they mostly emerge at radii < 100 AU. The fine-structure lines 
are only slightly affected with [C II] being slightly weaker and 
[C I] slightly stronger, owing to the greater charge exchange. We 
also tested the influence of the charge exchange with sulfur (C + 
+ S — > C + S + ), which had been traditional ly found to influence 
[C I] in PDRs (iHollenbach & Tielensll 19971 Section 4.4). For the 
dense PDR of the disk atmosphere including a density gradient, 
the influence of sulfur charge exchange is however found to be 
smaller than that for PAHs. 



Varying the PAH abundance in the inner 100 AU does not af- 
fect the atomic fine-structure lines. The high-/ lines of CO how- 
ever change by about a factor of two, if the PAH abundance is 
higher by a factor of five or set to zero. This is a result of either 
the higher or lower gas temperature caused by the contribution 
of photoelectric heating on PAHs to the thermal balance. 
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a.) Different outer PAH abundances (r > 100 AU) 
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b.) Different inner PAH abundances (r < 100 AU) 



x(PAH) - In: 
x(PAH) - In: 
x(PAH) - In: 



ISM, Out 
ISM. Out 
5 ISM. Out 



0.2 ISM 
0.2 ISM 
0.2 ISM 



X Observed Fluxes 




■A- 



x 
■ • 



10 



15 20 
CO - J llD 



2o 



30 



3 

3. 



a a 5 



« O. 



Fig. 12. Dependence of the integrated 
line fluxes on the PAH abundance, a.) 
The abundance of PAH outside 100 AU 
is varied, b.) The PAH abundance within 
100 AU is varied. 
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Fig. 13. Line fluxes of models with dif- 
ferent stellar radiation fields. 



4.4. The stellar input spectra 

The stellar UV radiation field is an important input parameter 
for the model, as both the heating and the chemistry depend 
on the input radiation field. We use observed and dereddened 
IUE/FUSE spectra. However, the dereddening is somewhat un- 
certain, particularly for shorter wavelengths. We thus also run 
models using a blackbody of 8000 K, 10500 K (T eff of HD 
100546), and 12000 K (FigureQ]). The blackbody was scaled to 
preserve the bolometric luminosity. Thus, the luminosity in the 



FUV band changed. While Lpuv/^boi ~ 0.27 for the observed 
FUV spectra, it is L F uv/£boi ~ 0.024 (0.095, 0.15) for the black- 
body of 8000 (10500, 12000) K. 

In Figure [13] the line fluxes of the models with the differ- 
ent input spectra are shown. Models with a lower FUV lumi- 
nosity Lfuv have less heating and thus the mid/high-/ CO and 
[OI] lines are weaker (Section O- The [CI]/[CII] ratio how- 
ever changes only slightly with the UV color. We note that both 
C and CO are photoionized/photodissociated at similar wave- 
lengths (A <1 100 A) and their rates thus depend similarly on the 
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UV color. For much colder UV spectra with very few CO dis- 
sociating photons, the gas temperature is lower owing to more 
efficient cooling by CO than for the atomic fine-structure lines. 
We conclude that for the lines discussed in this work, the UV 
spectrum affects the lines mostly by means of the photoelectric 
heating on PAHs/VSGs, which is less efficient for lower Lfuv- 
The [C I]/[C II] ratio remains however similar and the results re- 
lated to their flux are unaffected by an uncertainties in the dered- 
dening of the UV spectra. 

4.5. Dependence on the H 2 formation rate 

A major uncertainty in the chemistry is the H2 formation rate 
at high temperature. This rate may affect the line fluxes both in 
terms of the change in the H/H2 transition, which also affects the 
C + /C/CO transition and the thermal balance. The heating and 
cooling processes affected by H2 are cooling by H2 lines and 
heating by means of UV pumping followed by collisional deex- 
citation, formation, and photodissociation heating. Uncertainties 
enter the H2 formation rate in terms of the minimum grain size 
that determines the grain surface on which H2 may form, but also 
through the sticking coefficient and the formation efficiency. The 
sticking coefficient gives the probability that H atoms stick onto 
the surface rather than bounce back or evaporate before reacting 
to H2. This sticking coefficient an d the formation efficiency de - 
pend on various parameters (e.g. iHollenbach & Tielensl fl999). 
It may either decrease with the gas temperature when H2 for- 
mation on physisorbed sites dominates, or stay constant when 
at high gas temperature the neighboring chemisorbed H atoms 
react and leave the surface relatively bare. 

In this section, we calculate the representative model with a 
sticking coefficient times formation efficiency S X rj = 1, that 
is independent of temperature. We then compare the result to 
our "standard" H? f ormation rate using the S x rj reported by 
Cuppen et al. (201 0j), in o r der to capture the theoreti cal work of 
Cazaux & Tielensl (|2004|) . ICuppen & Herbstl (l2005h as well as 
the observational evidence from lHabart et al.l (120041) . Figure FBI 
shows the line fluxes calculated with the models assuming differ- 
ent formation rates. The inset in that figure gives the temperature 
dependence of the H2 formation rate. 

Neither the low-/ to mid-7 CO lines nor the atomic fine- 
structure lines depend on the H2 formation rate. The high-/ lines 
of CO, however, increase by a factor of ~ 3 when the sticking co- 
efficient is assumed to be independent of gas temperature. This is 
consistent with in the inner disk (r < 70 AU, z/r > 0.15) hydro- 
gen being fully molecular (see Figure [3] and HJ, which increases 
both the gas temperature and CO abundance in this region. We 
conclude that the predictions of the CO ladder are robust for the 
low-/ and mid-/ lines, although some differences due to the un- 
certain H2 formation rate are expected for the high-/ lines with 
/up > 25. 

5. Discussion 

5.1. Comparison with observations 

How do the results of the models compare to the observed 
fluxes? Studying the outer radius, gas/dust ratio, and carbon 
abundance in Section l4~Tl (Figure[8Tl. we find that the mid-//high- 
/ CO lines can be reproduced by either a higher gas/dust ratio 
and a lower 6c or a lower gas/dust ratio and a higher 5c- The 
low-/ CO and [O I] lines emerge from the outer part of the disk, 
hence scale with the outer radius of the disk. It is thus possible to 
reproduce these lines for both higher gas/dust and lower gas/dust 



ratios, if the outer radius is adjusted accordingly. For a gas/dust 
ratio of 20 and 6c = 0.5, a 400 AU disk better reproduces the 
[OI] fine-structure lines as the 200 AU disk which underpro- 
duces the flux. For a gas/dust ratio of 100 and <Sc = 0.05/0.1, a 
smaller disk better reproduces the [O I] observations, but a 400 
AU disk does not considerably overproduce the emission. 

For both gas/dust ratios and outer radii, the modeled low-/ 
lines agree reasonably well with the observations. Independent 
of the outer radius, the [C II] line is underpredicted for the mod- 
els that fit the CO ladder and [OI]. For gas/dust = 20, 6 C = 0.5, 
and r out = 200 AU, the [C II] is only underproduced by about a 
factor of three, but that model overpredicts the [C I] line by about 
a factor of four. However, for gas/dust =100 and 6c =0.1, both 
models with r out = 200 AU and 400 AU underproduce [CII], 
but are within the upper limit to [CI]. The same is true for the 
model with 6c = 0.05, which however more closely reproduces 
the CO ladder. Thus, none of the models are able to simultane- 
ously reproduce both the [C II] line and the upper limit to [C I] . 
We conclude that only models with a higher gas/dust ratio and 
low 6c ^ 0.1, are able to simultaneously reproduce the CO lad- 
der, [O I], and the upper limit to [C I] but underpredict [C II] by 
a factor depending on the outer radius of the disk. 

Changing the dust opacities from Ry = 3.1 to Ry = 5.5 does 
not change these conclusions (Section l4~2l Figure[9]). For models 
with 0.1-1.5 /mi grains and a much lower UV/visual extinction, 
the models with gas/dust =100 overproduce both the CO ladder 
and [OI]. For gas/dust = 10 and 6 C = 0.1, both the CO ladder 
and the [O I] lines are fit reasonably. The [C II] line is underpro- 
duced, but [CI] is within the upper limit for both radii. This is 
similar to our results for the models with a Ry =5.5 dust opacity 
law, gas/dust = 100, and 6 C = 0.05 /OA. 

Decreasing the PAH abundance in the outer disk (radii > 100 
AU) mostly affects the [CI]/[CII] ratio and a model without 
PAHs in this region starts to approximate both lines (Section l4~3l 
Figure [T2l. A higher PAH abundance in the inner disk (< 100 
AU) increases the flux of the mid-/ and high-/ CO lines and 
yields a slightly closer agreement with the high-/ lines, despite 
slightly overproducing some mid-/ lines. Given the large uncer- 
tainties in the gas temperature (Section 15.5b , we however con- 
clude that these changes are not significant. 

Using an H2 formation rate with a sticking coefficient that 
is independent of temperature slightly increases the high-/ lines 
of CO, which are in closer agreement than using our "standard" 
sticking coefficient that decreases above a certain temperature. 
However, as discussed before, the uncertain gas temperature 
most affects the high-/ lines and we thus do not consider the 
change in the high-/ lines with H2 formation rate as significant. 

We conclude that the CO ladder together with [O I] and the 
upper limit to [C I] can be reproduced by models with a high 
gas/dust ratio and low <5c- These models however underproduce 
[C II] . Models without PAHs in the outer disk, improve the fit to 
[C I] and [C II] but there is still a factor of two discrepancy for 
both lines. 

5.2. The [C\]/[C II] ratio and the gas carbon abundance 

Given the difficulties in simultaneously reproducing the [CII] 
line and the [C I] upper limit, it is interesting to study the depen- 
dence of the two lines on the parameters independent of the other 
lines. The modeled fluxes of the [C I] 370 /im line and the ratio to 
the [CI] 158 fim line are shown in Figure [15] which summarizes 
Figures[8][T0]and[l2l 

The dependence on r out , gas/dust, and 6c is shown in the up- 
per panel of Figure Q3] for models with the different dust opaci- 
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Fig. 14. Fluxes for models assuming dif- 
ferent sticking coefficients/formation ef- 
ficiencies to calculate the H2 formation 
rates. The temperature dependence of 
the rate is given in the inset. 
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b.) Different PAH abundances 
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Fig. 15. Modeled flux of the [C I] 
370 fim line and the line ratio 
of [C I] 370 fim / [C II] 158 fim. 
Note that in panel a.), gas/dust 
= 10 is given for the 0.1-1.5 fan 
opacities, and gas/dust = 20 for 
i?v = 3.1,5.5. The symbol 6c is 
the fraction of carbon in the gas- 
phase and r out is the outer radius 
of the disk. In panel b.), x(PAH) 
is the PAH abundance relative to 
the ISM abundance within and 
outside 100 AU. 



ties used in this work. The [C I] 370 fim line and the [C I] 370 fim 
I [C II] 158 fim ratio depend similarly on r out , gas/dust and 5c for 
the different opacity laws. Exceptions are the large grain models 
owing to the UV photons penetrating further into the disk, thus 
dissociating precursors of CO, such as CH and CH2, and affect- 
ing C/C + . The [C I] 370 fim line varies by about a factor of four 
between the models with outer radii of 200 AU and 400 AU and a 
similar factor between 6 C = 0.1 and 0.5. The [CI] 370 fim /[C II] 
158 fim line ratios are almost independent of both the outer ra- 
dius and carbon abundance, but depend slightly on the gas/dust 
ratio, again since CH and CH2 are more likely to be converted 
into CO than photodissociated at higher density, thus decreas- 
ing the [C I] fluxes. Thus, by increasing the gas mass by a factor 
of five from gas/dust of 20 to 100 only increases the [CI] by 
about a factor of two. This reduces the [C I] 370 fim /[C II] 158 
fim line ratio by about a factor of two. For the larger grains, this 
effect becomes much stronger because the lower UV extinction 
moves the C + /C transition deeper into the disk, to higher den- 



sities. We note however, that the models with large grains and 
gas/dust = 100 overproduce CO and [OI] (Sectionl4~2l). 

The PAH abundance within 100 AU does not affect the [C I] 
and [CII] fluxes. However, if PAHs are absent in the outer part 
of the disk (> 100 AU), the [CI] fluxes decrease by about a 
factor of two, and the [CI] 370 fim /[CII] 158 fim line ratio 
decreases by about a factor of three. As discussed before, this 
is due to the lower charge exchange from C + to C by PAHs. 
Similarly, increasing the the PAH abundance from 20 % ISM to 
ISM abundances leads to a stronger charge exchange and thus a 
larger line ratio. 

In spite of all these parameter changes, the observed [C II] 
158 fim line and the observed line ratio [CI] 370 fim / [CII] 
158 fim cannot be reproduced by our models. Our closest match 
is found for the model without PAHs in the outer part of the 
disk (Figure [T2l. but that model still differs by a factor of about 
two for both lines. The [CI] 370 fim / [CII] 158 fim line ratio 
is surprisingly independent of the parameters, which can be un- 
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derstood by a similar region of origin of the two lines found in 
Section [141 

5.3. The origin of the [C \\] emission? 

Does the observed [CII] 158 /mi emission really originate from 
the disk? Since C + is the main state of carbon in diffuse clouds 
and owing to the low critical density of the 158 /mi line, how 
much could a tenuous remnant envelope or a diffuse foreground 
cloud contribute to the emission? We note that the discussion 
here assumes that the emission is centered on the source, within 
a ~ 10" (1000 AU) diameter region, since extended emission 
has already been subtracted (Section [Z2b . 

If the extinction towards the star results from material only 
seen in the central pixel, how much could that contribute to the 
[C II] emission? The maximum amount of emission to be ex- 
pected from the foreground can be calcul ated from the fore - 
ground extinction of Ay = 0.36 mag dArdila et alJ 120071) . 
Following Appendix lB.il we obtained a flux of ~ 2.5 x 10~ 16 
W irT 2 in the Herschel beam , using a conversion of 1.9 x 10 21 
cm 2 A v l (Bo hlin etalJfl978h . a carbon fractional abundance of 
about 1.3 x 10 -4 , and a gas temperature in the foreground above 
the upper level energy of the [CII] 158 /mi line (91 K). This 
would already exceed the flux reported in Table [T]by 60 %. 

How much emission could come from a compact diffuse 
remnant env elope? Scattered lig ht is detected out to a distance 
of 1000 AU dArdila et alJl2007h and could belong to either the 
disk or a remnant envelope. If a tenuous remnant envelope is ex- 
posed sufficiently to UV radiation, most of its carbon is in the 
form of C + . This nebula cannot contain CO, since this would fill 
in the central emi ssion of the observed double-peak line shape 
dPanic et al.ll2010h . While the density structure of such a rem- 
nant envelope is not known from observations, assuming a ring 
size in the range ~ 400- 1000 AU at a density of ~ 10 4 cm~ 3 and 
a thickness of ~ 800 AU would yield a flux of up to 1.4 x 10~ 17 
W irT 2 . Thus, both a foreground cloud and a remnant envelope 
might well contribute considerably to the [CII] emission. We 
note that no [CI] emission is expected from this component, 
since for typical columns of C towards diffuse clouds (< 3 x 10 15 
cirT 2 . lJenkins & Shavall979l) . the [C I] 370 fjm line is below the 
detection limit of APEX. 

This rough estimate shows that both foreground and remnant 
envelope could add substantially to the [C II] emission, making 
the association of the observed flux with the disk not obvious. 
This question can however only be addressed with velocity re- 
solved spectra of [CII] 158 /mi, obtained with the HIFI instru- 
ment onboard Herschel. 

5.4. The carbon budget of the disk atmosphere 

Our results indicate that for normal gas/dust ratios, the car- 
bon abundance is low at 5q 5= 0.1. How does the gas-phase car- 
bon abundance change as the gas evolves from diffuse clouds to 
dense clouds, collapsing envelopes, and eventually disks? At the 
earliest stage, in diffuse clouds, the carbon gas-phase abundance 
can be measured directly by t he UV absorption of C II, without 
referring to a dust model (e.g. lSofia et al.l2004l.l201 1|) . The mea- 
sured amount of volatile carbon, which is not bound in dust, is 
constrained to be between ~ 40 % and 60 % (6 C = 0.4 to 0.6), 
with some differences between the measurements of the strong 
1334 A feature dSofia et al.ll20 lib and the weak intersystem fea- 
ture at 2325 A dSofia et al.ll2004l) . Assuming that the refractory 



material does not change as the gas evolves to cores and finally 
disks, these micron-sized carbon grains quickly settle down to 
the midplane of the disk. Thus, if a process removes gas from 
the upper disk, that population of carbon dust is not removed 
from the disk. 

The volatile carbon however undergoes significant changes 
in the course of this evolution: from diffuse to dense clouds, 
the main reservoir becomes gaseous CO rather than C + . In 
pre-stellar cores, it be comes CO ice rather than CO gas 
(lBergin&Tafaliall2007l) . From pre-stellar cores to disks, it ei- 
ther stays in the form of CO ice, adsorbs or desorbs from grains 
multiple tim es, or is transformed into o ther species such as CO2 
or CH 3 OH dVisser et al.ll2009bl 1201 lh . If gas is removed from 
the disk by some process, leading to a smaller gas/dust ratio, the 
volatile carbon is likely to be removed at a rate that is propor- 
tional to the quantity of gas. Thus, within the gas, the amount 
of volatile carbon should remain the same. Inferring a value of 
6c that is much lower than 0.4 to 0.6 as we find here means that 
a fraction of volatile carbon has been turned into non-volatile 
material in the course of the evolution from diffuse clouds to 
disks. Towards hot c ores, where CO ice is supposed to evapo- 
rate fr om dust grains. lAlonso-Albi et alJ d2010l) and lYildiz et al.l 
d2010h found evidence of low CO hot core abundances, suggest- 
ing some transformation of CO on the surfaces of the grains into 
other volatile or non-volatile carbon species. For the warm disk 
discussed here, CO ice cannot be the main carbon reservoir. In 
addition, no strong CH3OH emission is in general observed in 
disks, ruling out this possible carbon reservoir. However, some 
"mild" ph otochemistry could lead to more complex, less volatile 
organics dOberg et al.ll20 09a). Some of this less volatile organic 
material could become simi lar to the so-ca lled "CHON parti- 
cles" detected in comets dKissel et al.l Il986h after further pro- 
cessing in the disk. 

The upper layers of the disk atmosphere, probed in this work, 
more likely have gas/dust > 100 rather than < 100, since settling 
increases the gas/dust ratio in that region. In addition, as we keep 
the dust density structure in our work fixed, a small gas/dust ratio 
would mean that the relative fraction of carbon locked into dust 
( VSGs/PAHs) increases for small gas/dust ratios, which actually 
implies that the total carbon abundance is higher than the solar 
abundance. These arguments drive us to the conclusion that the 
warm gas in the HD 100546 atmosphere more likely consists of 
a large gas/dust ratio with a considerably smaller abundance of 
volatile carbon than in the diffuse ISM. 

5.5. Uncertainties in r gas 

An important caveat of the models presented here is the un- 
certainty in the derived gas temperature. On the surface of the 
disk where the gas and dust temperature decouple, the equilib- 
rium of different heating and cooling processes govern the gas 
temperature. As noted before, the thermal balance is a coupled 
problem together with the chemistry as line cooling (in [CII], 
[OI], CO, ...) depends on the abundances. From PDR mod- 
eling, the thermal balance problem is known to be very sensi- 
tive to the assumptions made, such as the method used for the 
line radi ative transfe r or th e H2 physics and chemistry imple- 
mented. Rolli g et alJ d2007l) compared the gas temperature de- 
rived by different PDR models assuming the same chemical net- 
work. They found a relatively large scatter, in particular for the 
model with a high density (« = 10 s 5 cm 4 ) and a strong UV 
field {x — 10 5 ) corresponding approximately to the upper atmo- 
sphere of the outer disk. The scatter between the models is about 
a factor of three in r gas (Figure IA.7I ) and similar for the calcu- 
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lated fine-structure line fluxes. In our present study, we found 
the high-/ lines of CO to be more sensitive to the temperature 
(Section [3~3l and we thus calculated the CO ladder from the CO 
abundance and gas temperature obtained from the different PDR 
codes. The scatter in those line fluxes was found to be about 1.5 
dex when outliers were removed. The chemical network, which 
has been taken to b e the same th an the one used in the bench- 
mark study bv l Rollig et al.l d2007l) . also introduces uncertainties. 
Vasv unin et al.l d2008l) performed a Monte Carlo sensitivity anal- 
ysis and determined the uncertainty in the CO and C + column 
density through a disk to be a factor of two or less. 

For our work, we thus consider an agreement of a factor 
of two between observations and model as a good agreement. 
Given that many model parameters remain free even for this 
well-studied disk, the observed fluxes can be reproduced to a 
level that is much better than a factor of two. We however re- 
frain from trying to achieve so, as we feel that given the large 
uncertainties, these results would not be trustworthy and provide 
a misleading picture. 

Given the considerable uncertainty and disagreement be- 
tween the models, which of our conclusions is robust? From 
the considerably higher gas temperature needed to reproduce the 
high-/ CO lines, together with all PDR models agreeing that the 
gas temperature indeed decouples from the dust temperature at 
these high densities and UV fluxes, we conclude that there must 
be a warm gas atmosphere. We also consider the conclusions 
made in Section|4]as robust, since they mainly represent relative 
changes and trends that are more reliable than absolute fluxes, 
as long as the correct main heating and cooling processes domi- 
nate the thermal balance. Similarly, finding trends in the depen- 
dences of the line fluxes on model parameters, which can then 
be observationally tested using a larger sample of disks, is more 
promising than fitting the absolute fluxes. 

5.6. The shape of the CO ladder 

The high-/ CO lines emerge from the upper, inner atmosphere 
of the disk, from regions where heating by FUV radiation leads 
to a gas temperature that is much higher than the dust temper- 
ature (Section I3.4t . However, what does the shape of the CO 
ladder and in particular the mid-/ and high-/ CO lines tell us? 
In Figure [16] we show the rotation diagram of the CO ladder 
fro m both observations and th e representative model (Section 
0. Golds mith & Langerl (1 19991) provide a detailed introduction 
to the rotation diagram method. The Y axis of the diagram is de- 
fined to include the solid angle of the emission, because the CO 
ladder is thought to emerge from different regions (Appendix 

E3 . 

The observed CO ladder can be represented by three tem- 
perature components: A cold component with a rotation tem- 
perature of T rot ~ 60 K (/ up = 3 - 7), a warm component 
of T m ~ 300 K (/ up = 14 - 25), and a hot component with 
/rot ~ 750 K (/ up = 25 — 30). We did not include the upper lim- 
its to the / = 23 - 22 and / = 31 - 30 lines in the calculation of 

/rot- 

The low-/ lines are optically thick. This leads to an increase 
in the flux compared to the optically thin conditions (Appendix 
IB .2) . Owing to the opacity effects, the temperature obtained 
from the cold component is thus not necessarily a kinetic tem- 
perature even though the lines are in thermal equilibrium. The 
model predicts that lines up to E u ~ 750, / = 16 - 15 are op- 
tically thick (Section [5). The modeled CO ladder indeed shows 
a curved shape, up to about this upper level energy, which is 
characteristic of optical depth effects. However, a larger emit- 



ting area, leading to a larger dQ. can also increase the line flux 
in a similar way. The emitting region dQ. contributing more than 
50 % to the total emission (obtained from the 75 % radii given 
in Table [2} increases from / = 16-15to/ = 3- 2bya factor 
of about eight, corresponding to a shift of AY ~ 2 in the rotation 
diagram. 

The optically thin mid-/ and high-/ lines obtained from the 
model, can be represented well by a single temperature com- 
ponent of 393 K, which is somewhat higher than the observed 
value. The difference is however small in comparison to the 
large uncertainty in the gas temperature. The hot component of 
r, ot ~ 750 K is not reproduced well by the representative model. 
A higher H2 formation rate, leading to a higher H2 abundance 
in the hot, inner gas, thus greater heating and more CO, yields 
an additional hot component. This model more close reproduces 
the observations of the hot component (Section l4~5l >. 

How is the CO ladder affected by the molecular excitation? 
To test this, we compared the CO ladder from the representa- 
tive model with the same model, but assumed that the CO level 
population is in local thermal equilibrium (LTE). While the op- 
tically thick low-/ to mid-/ lines are very similar, the optically 
thin lines with higher / are stronger by a factor of about nine 
(AY ~ 2) in flux. The critical density of these lines is on the 
order of up to a few times 10 6 cirT 3 (Table [2}. In the region of 
the main emission ( Section [3~4l i. the density is similar or higher 
than the critical density. The calculated level population in this 
region is similar to the LTE population, to within ~ 10 % for 
/ = 16 — 15 and to within a factor of two for / = 30 - 29. The 
increase in the fluxes when assuming LTE is explained by the 
tenuous and hot gas at larger radii also contributing to the emis- 
sion, as in the case of for example the "warm" CO finger in the 
outer, upper disk, discussed in Section [3~!2l 




500 1000 1500 2000 2500 



E up (K) 

Fig. 16. Rotation diagram of the modeled and observed CO lad- 
der. 

We conclude that different effects govern the shape of the 
CO ladder including: ( i) The temperature of the emitting gas and 
the density, (ii) Optical depth effects at low-/ to mid-/ lines. 
( Hi) For optically thin lines, the number of emitting molecules 
(~ N(CO)dQ) and for optically thick lines, the emitting area 
(~ dQ). Depending on the transition, different effects dominate. 
In addition, different temperature, column density, and density 
structures might reproduce the same CO ladder. 
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6. Conclusions and outlook 

We have used a new detailed physical-chemical model to calcu- 
late the atomic fine-structure emission ([CI], [CII], and [OI]) 
and CO lines from CO J — 3 - 2 to J — 30 - 29 from a proto- 
planetary disk around a young Herbig Ae star. We have applied 
the model to the well-studied HD 100546 disk, but the results 
are more generally applicable. We have compared our results to 
recent observations of the object and studied the dependence of 
the line flux on different parameters. Our main results are: 

- The high-/ emission of CO can only be reproduced when 
the gas temperature in the inner disk is considerably higher 
than the dust temperature. A model with the gas temperature 
set to the dust temperature underpredicts the high-/ lines by 
orders of magnitude (Section l3~3l . The UV radiation is able 
to heat the gas sufficiently (Section lXTT i. 

- Both models with a gas/dust =100 and 6c = 0.05 (volatile 
fraction of carbon 6c) and models with gas/dust = 20 and 
6c = 0.5 are able to reproduce the CO ladder and the [OI] 
lines. Models with a high gas/dust ratio overpredict the upper 
limit of [C I] but approximately fit [C II]. Models with a low 
gas/dust ratio, however agree with the upper limit of [C I] but 
underpredict [CII] (Sectiongj). 

- None of the models reproduce [CI]/[CII], but models with 
higher gas/dust ratio have a lower ratio of these lines, closer 
to the observed value. The line ratio is also smaller in models 
with larger grains that allow UV radiation to penetrate deeper 
into the disk (Section l4~2ll5.2l i and models without PAHs in 
the outer disk, at radii > 100 AU. The ratio depends less on 
the stellar radiation field (Section l4~4l >. 

- A substantial fraction of the [CII] emission may emerge 
from a remnant envelope or a compact foreground (Section 
15.31 1. hence we prefer a solution with gas/dust ratio of 100, 
but a small fraction of volatile carbon (Section I5.41 i. This 
means that a considerable part of the carbon is locked into 
more refractory carbon-bearing species and that the carbon 
partitioning between gas and dust has evolved from diffuse 
clouds to protoplanetary disk. 

- The highest-/ emission detected towards HD 100546, CO 
/ = 30 - 29, emerges from radii of ~ 20 - 50 AU. Mid-/ 
lines, e.g. CO / = 16 - 15, are emitted at radii of ~ 40 - 90 
AU and low-/ lines trace the outer disk (Section [3.41 . 

- The high-/ lines of CO (/ up > 14) are predicted to be much 
broader than the low-/ lines, which are observable from the 
ground (Section [3~5l l. 

- Changing the dust opacity law to very high values of Ry 
allows UV radiation to penetrate considerably deeper into 
the disk. The resulting line fluxes for a fixed gas/dust ra- 
tio are much larger. For 0.1-1.5 /mi grains, gas/dust = 10 
with 6c = 0. 1 reproduces the observations similarly well as 
gas/dust = 100 and 6c = 0.05, if a Ry = 5.5 dust opacity law 
is used (Section l4~2l . 

- Assuming a temperature-independent sticking coefficient for 
the H2 formation only affects the high-/ lines, changing the 
line fluxes by a factor of a few. The uncertain H2 forma- 
tion rate is thus not a critical parameter for the CO emission, 
except for the highest-/ lines detected towards HD 100546 
(Section Pl E3 . 

- The shape of the CO ladder in the rotation diagram is gov- 
erned by the combination of excitation effects, optical depth, 
and different emitting regions of the CO lines (Section [5.61 . 

With this paper, we have introduced a new physical-chemical 
model and applied it to a protoplanetary disk. The model is based 



on our previous work dBruderer et alJl2009al 1201 Oh and calcu- 
lates the gas-temperature structure in a self-consistent way with 
the chemistry and the excitation of the molecules. The model has 
a wide range of applications and we plan to use it in the future 
to investigate both T Tauri disks as well as the outflow walls of 
the envelopes of YSOs, and to compare the results with Herschel 
observations. The model also can be relatively easily adopted to 
calculate three-dimensional structures and is well-suited to the 
analysis of upcoming high-resolution ALMA data. 
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Appendix A: Details of the models and benchmark 
tests 

The purpose of this appendix is twofold. First, it presents de- 
tails of the implementation of the model and second, it com- 
pares of results with other c odes and benchmark pr o blem s. The 
model is based on those of lBruderer etal.1 (l2009blH 120101) and 
lBrudererld201Q ). although it contains several improvements. The 
structure of the model is summarized in Fig. IA.ll The modeling 
process starts with a gas and density distribution e.g. obtained 
from an analytical prescription. Different modules then calcu- 
late the dust temperature and UV radiation (Sect. IA.2K abun- 
dances of atomic or molecular species (Sect. IA.3l l. and the level 
population of molecules or atoms acting as coolants (Sect. lA~4l i. 
The gas-temperature is then calculated from the balance between 
heating and cooling processes (Sect. lA3T l. As the abundance and 
level population enter the heating and cooling rates, abundances 
and level populations need to be re-calculated until convergence. 
Finally, spectral lines can be calculated and compared with ob- 
servations. 

To carry out the calculation, the modeled region is divided 
into individual cells. In this present study, we assume a axisym- 
metric structure. Since all geometric information is contained in 
one separate module, the code can in principle also be run for 
spherically symmetric or three-dimensional structures. 
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Fig. A.l. Modeling flowchart 



The current study applies the model to a protoplanetary disk, 
but the model can also be used in the context of protostellar en- 
velopes, for example to calculate the molecular emission from 



r, 



UV irradiated outflow walls. For future applications, we thus 
present all modules in this appendix, even though they are not 
used in the current work (such as the calculation of the dust tem- 
perature). 

A.1. Dust radiative transfer 
AAA. Dust temperature 

The dust temperature in steady state is obtained from the equi- 
librium between energy absorbed and emitted by dust 

bs = K A I A dAd£l = 4-n J K A B A (T dusl )dA = r em i t , (A.l) 

where the dust opacity per unit mass is k a , the incident radia- 
tion field is I A , and the Planck function for the dust temperature 
is fi,i(7dust)- The radiation field consists of contributions from 
both the stellar radiation field and the ambient radiation by dust 
emission. Thus, different positions in the model are coupled by 
means of the radiative transfer equation. We neglect the "back- 
warming" of dust by hot gas and further assume a unique dust 
temperature for different sizes of dust. 

In this present study, we implement the Monte Carlo tech- 
nique proposed by Bior kman & Woodl (l200lh : photon packages 
with an energy distribution following the stellar radiation field 
are propagated until interaction (scattering or absorption fol- 
lowed by re-emission). Scattering only changes the direction 
of a photon package. To sample the scattering angle, we use a 
Henyey-Greenstein function, which defines the scattering angle 
by only one parameter, g = (cos(f?)). Absorption increases the 
amount of energy deposited into a cell. For a cell of dust mass 
M, the dust temperature T nevl after absorption of Nj photon pack- 
ages is given by 



Nj5L 
47rM' 



(A.2) 



with the energy per photon package 6L = L/N, where the stellar 
luminosity L is divided by the total number of photon packages 
N. When an absorption event increases the local temperature by 
Ar = - r bu St , the photon-p ackage is re-emitted at a wave- 
length given by (iBaes et al.ll2005l) 



n j ,xi K v B v (T Dust + AT)dv K v B v (T Dust + AT)dv 
Pydv = (Nt + 1) — r. Ni — r. 

J KyBy(T Dust )dV J KyBy(T Dast )dV 



.(A3) 



This probability distribution function (PDF) accurately cor- 
rects for the incorrect spectrum of previously reemitted pho- 
tons (at a too low temperature). However, it can sometimes be- 
come negative. In such c ases, the approximative PDF given in 
iBjorkman & W ood (2001) is used. This procedure is generally 
very efficient as it does not require any global iteration and the 
dust temperature increases to the correct value. In regions of very 
high optical depth, however, photons can get stuck owing to a 
very high number of absorption and re-emission proc esses. We 
thus im plement t he diffusion appro ximation derived in Robitaille 
(2010), based on lMin et all (|2009). Regions of very low optical 
depth can suffer too little absorptio n resulting in a noisy dust 
temperature (see e.g. Bruderer 2010, Fig. D.2). We implement a 
method proposed by Pint e et all d2006l) and measure the a mbient 
radiation field using the estimator given in iLucvl (119991) during 
the propagation of the photons. After all photons have propa- 
gated, the dust temperature can be improved using this value. 

To benchmark the dust temper ature calculation, we use the 
benchmark problems proposed in llvezic etal I( fl997h and com- 
pare our results with calculations obtained with the DUSTY 
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code dlvezic & Elitzurlfl9 97). The benchmark problem consists 
of a spherical cloud with two different density gradients (either 
oc r~ 2 or constant) and four different total radial optical depth 
at 1 fim (ti m = 1, 10, 100, or 1000). The input spectra is a 
black body and the dust properties follow a simple analytical 
law. Scattering is assumed to be isotropic. The agreement be- 
tween the codes is perfect, being closer than about a percent. 



1000 



r,(](j 



200 




10 100 
Radius (y = r/r Dust ) 



11 



1000 



r.oo 



t>-) p = poy 2 




10 100 
Radius (y = r/r Dust ) 



1000 



Fig. A.2. Results of a spherical benchm ark test for the dust tem- 
perature calculation dlvezic et al.l 1 1997b . The dust temperature 
obtained from this work is given by a red line, while the results 
from DUSTY are given by black signs. The solid line corre- 
sponds to T\ m , and the dashed, dotted, and dash-dotted lines 
to T\ = 10, 100, and 1000, respectively, a.) constant density, b.) 
density oc r~ 2 . 

As a second benchmark problem, we calculate the dust tem- 
peratu re in a protoplan etary disk with very high total optical 
depth dPinte et alJ200 9i). We calculate their benchmark problem 
for isotropic scattering with a midplane optical depth at 0.86 fim 
of r = 10 3 and 10 5 (their Fig. 9). In Fig. IA.3I we show the mid- 
plane temperature for the two models compared to the results ob- 
tained with MCFOST. The agreement is very good for the model 
with t = 10 3 . For the run with t = 10 5 , the agreement is again 
very good with deviations smaller than a few percent, even in 



the region shadowed by the dense inner rim. These deviations 
are comparable to those between the other codes participating in 
that benchmark study. 



a.) Midplane Temperature for r = 10' 




10 -4 

b.) Midplane Temperature for r = 10 




10 _1 10° 
R - R mm (AU) 



Fig. A.3. Dust temperatu re at midplane for the disk benchmark 
problem (Pinte et al. 20091). Th e result from our code is given in 
red, the result from MCFOST (Pin te et al.ll2006l) in black. 



A.2. UV field 

Since the calculation of the dust temperature does not explicitly 
compute the radiation field as a function of wavelength, the in- 
tensity at UV and optical wavelengths is calculated in a second 
step. This intensity is used for example to calculate photodissoci- 
ation rates. The calculation of the intensity uses the same method 
for the propagation of the photon packages as the calculation of 
the dust temperature. The radiation field is then calc ulated with 
the method proposed by Ivan ZadelhoffetaD d2003h . For each 
wavelength bin, N photon packages are propagated. All photon 
packages are assigned with an initial intensity Ia(0) = FaS IN, 
where Fa is the input flux and S the surface it passes through. 
This intensity then drops to Ia(s + As) = Ix(s)e l , when the 
package advances by As with an optical depth At,}. The mean 
intensity is calculated from 



1 v- 1 

J A = > lAiASj — 



At, 



(A.4) 



with the volume of the cell V and the sum taken over all photon 
packages passing through a cell. The distance and optical depth 
as an individual photon package i passes through a cell are given 
by Asj and At,! ,-. Since we calculate the mean intensity for a few 
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wavelengths only, the input intensity is obtained by averaging 
over the input (stellar) spectra in each indivi dual spectral band. 
Photons from the interstellar radiation field (lDraineiri978f) are 
also accounted for. 

As a benchmark test, we calculate the mean intensity in a 
spherical cloud, switching scattering off. We compare the cal- 
culated intensity to the analytically obtained intensity and find 
agreement to better than one percent. The implementation of the 
scattering follows the calculation of the dust temperature and 
was tested in the previous paragraph. 

A.3. Chemistry 

Abundances of molecular and atomic species are obtained from 
the solution of the rate equations 

dn(i, t) \-i 

j fl 

with the abundance n(i, t) (cm -3 ) of a species i at time t. The rate 
of destruction and formation of a species is given by the coeffi- 
cients kjj and kiji, respectively. In addition, three-body reactions 
can be entered, but our current network does not contain any. The 
set of coupled, stiff, and non- linear differential e quations I A. 5l is 
solved by the DVODE solver dBrown et alJ[l 989) for the tempo- 
ral evolution of the species. We verify that the solution complies 
with the conservation of elements and charge. We note that the 
DVODE solver is wi dely used for applications in astrochemistry 
and well-tested (e.g. Neiad 2005). To obtain steady-state abun- 
dances, Eq. IA.5l is solved for dn(i, t)/dt = using a globally con- 
vergent Newton-Raphson method. To close the system of equa- 
tions Eq. IA.5I some rate equations are replaced by conservation 
equations. If the steady-state solver fails to converge, the abun- 
dances are obtained from the time-dependent solver running to a 
high chemical age (> 10 Mio. years). 

To test numerical aspects of the solver, we calculate abun- 
dances in an isothermal slab of constant density with UV irra- 
diation. A small chemical n etwork from the PDR comparison 
study by Rolli g et al.l d2007l) is used for this test. We first cal- 
culate their problem F4 (for a density of 10 5 ' 5 cm -3 , UV field 
of x — 10 5 , and a temperature of 50 K) with the steady-state 
solver. In Fig. IA.4I we give the abundances of H, H2, O, O2, 
C + , CO, and electrons and compare them to the results obtained 
from PDR codes. Our code to within their scatters agrees with 
the other codes. The scatter in the H2 abundance is due to a dif- 
ferent implementation of some reactions, such as H2 dissocia- 
tion including self-shielding. We also calculated the abundances 
with the time-dependent solver running to an old chemical age 
and found an agreement to better than a factor of 10~ 4 with the 
steady-state solution. 

A.3.1. Chemical network 

The chemical network used in our work contains 10 elements, 
109 species, a nd 1492 reactions. W e us e the same e l ementa l 
composition as lJonkheid et al.l d2006h and lWoifke et alJ (|2009a), 
which is summarized in Table I A. fl The species contained in the 
network are given in Table lA.2l We note that W 2 is vibrationally 
excited H 2 . PAH , PAH + , and PAH" are neutral, positively, and 
negatively charged PAHs respectively and PAH:H denotes hy- 
drogenated PAHs. Species frozen-out onto dust grains are shown 
by JX, for example JCO is CO frozen-out. 

The chem i cal re a ction network is ba sed on the studies 
bv lDotv et ail d2004l) . IStauber et all d2005l) . and iBruderer et alJ 



Table A.l. Elemental composition. a(b) means a x 10*. 



Element 


number fraction 


H 


1 


He 


7.59(-2) 


C 


2.70(-4)<5 c a 


N 


2.14(-5) 


O 


2.88(-4) 


Mg 


4.17(-6) 


Si 


7.94(-6) 


S 


1.91 (-6) 


Fe 


4.27(-6) 



Notes. (n) depending on the amount of volatile carbon, see main text. 

Table A.2. Species contained in the chemical network. See 
Section lA.3.1l for an explanation of the species. 



(2009b), but contains some extensions. The following reaction 
types are included: 

1. ) H 2 formation on dust 

Fo r the formation of H2 o n dust grains, we implement the rate 
bv lSternberg & Dalgarnol dl989l) . scaled t o the appropriate grai n 
size. We use the sticking coefficients of Cuppe net alJ ([2010). 
Since the gas/dust ratio can vary over a model, rates involving 
small dust are scaled to the local gas/dust ratio. 

2. ) Freeze-out, thermal, and non-thermal desorption 

Freeze-out and the rmal o r non- thermal desorption are imple- 
mented following IVisserl (l2009h and references therein. The 
thermal desorption rates account for the abundance of the 
ice on the surface and change from zeroth to first order be- 
havior when the_jce contains less than a monolayer of ice. 
As in Ivisseil 42009), we use the same pre-ex ponential fac- 
tor for all species exc ept for H2O dFraser et alj|2001|) and CO 
(Biss chop et alJ 120061). Binding energies f ollow lAikawa et al.l 
( 1 19971) and lSandford & Allamandolaldl993h . 

Non-thermal evaporation by UV photons or cosmic -rays de- 
pends on a yield indicating the number of molecules desorbed 
per grai n per inciden t UV photo n. We assume the same y ields 
used bv lVisseri (120091) . following lOberg et all d2007ll2009bl) . To 
simulate cosmic-ray induced desorption, a small background UV 
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Fig. A.4. Result of the chemistry benchmark for a slab with UV irradiation (Rollig et al. 2007, problem F4). The result of our code 
is given by the red line, the results of the PDR codes by gray lines. 



flux corresponding to 10 4 photons cm 2 s 1 (Shen et al . 2004) is 
added to the stellar UV flux. 

3.) Hydrogenation of simple species on ices 

In addition to freeze-out and evaporation, we implement the 
grain-surfac e hydrogenat ion of simple species. We f ollow the 
approach of lVisseri d2009l) and lHoUenbach et alj d2009h . The hy- 
drogenation leads to saturation of species frozen-out, from JC 
-> JCH -> JCH 2 -> JCH 3 -> JCH4, JN -> JNH -> JNH 2 -> 
JNH3and JO — > JOH — > JH2O. Other grain-surface reactions to 
build up more complex species are not accounted for, but the 
focus of our work is on simple species. 

5. ) Gas-phase reactions 

Most reactions in our network (1 132 out of 1492) are pure gas- 
phase reactions, such as ion-neutral, neutral-neutral, or charge 
exchange reactions. Our n etwork contains all r eactions from 
the UMIST 2006 database dWoodall et al.ll2007t) involving the 
species given in Ta ble IA.2I Details of the implementation are 
given in feruderer et all (l2009bT) . 

6. ) Photodissociation 

Photodissociation rates are calculated from the local radiation 
field using the cross-sections for discrete and continuous absorp- 
tion given by [v an Dishoeck ( 1988), van Di shoeck et alj (|2006), 
and lvan Hemert & van Dishoeck! (I2008fl For the dissociation of 
H2, CO, and the ionization of C, we include the effects of self- 

3 www.strw.leidenuniv.nl/~ewine/photo/ 



shielding by reducing the unshield ed rate with a self-shield ing 
factor. For H 2 , the factor given in iDraine &Bertoldild 19961) is 
used. F or CO, we implement the factors given in IVisser et alj 
d2009al) and for C, the factor of lKamp & Bertoldil (120001) is used. 
We however limited the shielding of C by H 2 to 0.5, owing to the 
overlap of H 2 lines with the wavelengths at which carbon is pho- 
toionized. 

The self-shielding rates depend on the column densities of 
H 2 , C, and CO towards the UV source. We calculate these col- 
umn densities together with the UV field (Sect. IA.21 i by averag- 
ing over the densities of all photon packages passing through a 
cell, weighted by their intensity. 

7. ) X-ray induced processes 

As demonstrated by iBruderer et al. (2009b}, the influence of 
X-rays on the chemistry can be accurately approximated by 
secondary ionizations induced by fast photo-electrons from H 2 
and H ionization. We follow their approach and calculate the 
X-ray ionization rate using the cross-sections given in their 
Appe ndix A. The rates for the secondary processes are taken 
from Stau ber et al.1 d2005l) . We note that the cosmic-ray-induced 
photodiss ociation rates are sca led-up to the X-ray ionization rate 
following lsFauber et al.1 d2005l) . 

8. ) Cosmic-ray-induced reactions 

Rates for direct cosmic -ray-induced ionization processes are 
taken from the UMIST 2006 network. Cosmic-ray-induced pho- 
toionization rates are also taken from that databa se, with the ex - 
ception of CO, H 2 , and H, where we follow Stauber et al. (2005). 
For the CO dissociation, we use the rates of lGredel et al.l (11987) 
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imple mented with the fitting equation provided in lMalonev et al.l 
( 1996]). We account for the ionization of C O, H 2 , and H by the 
decay of excited He (2'P state, 19.8 eV, see lYanlll 9971) . 

9.) PAH/small grain charge exchange/hydrogenation 

We assume that PAHs and small grains are well-mixed with the 
gas and thus not scaled to the gas/dust ratio. For the photoion- 
ization, c harge-exchan g e, and recombination, we use the rates 
given in Wolf ire et al.l (120031) with a collision rate parameter 
0pah = 0.5. For species heavier than hydrogen, the recombina- 
tion rates are scaled to 1 / yjm. dmn , where m amu is the mass of the 
molecule. 

The for mation of CH + and H 2 on PAHs is implemented using 
the rates of iJonkheid et al.1 (120061) . We include the reactions 

i0 



H + PAH U -> PAH : H 
H + + PAH -> PAH : H 
H + PAH : H -> PAH + H 2 
C + + PAH : H -> PAH + CH + 



(A.6) 
(A.7) 
(A.8) 
(A.9) 



Reactions IA.61IA.7I and lA.8l run at 10 9 cm 3 s while lA.9l runs 
at the same speed as the recombination of C + with PAH . 

10.) Reactions with H* (excited state H 2 ) 

Vibrationally excited H 2 is included into the chemistry with 
a two-level approximation. H* 2 denotes H 2 in a vibrationally 
excit ed pseudo-lev el with an energy corresponding to 30163 
K dLondonl [l978h . The energy stored in the vibrational 
state can be used to over come an a ctivation barrier (e.g. 
ISternberg & Dalgarnolll995l) . Following iTielens & Hollenbach 
(119851) . the rate coefficients of a species with Hi can be approx- 
imated using the rate for H 2 by reducing the exponential factor 
by 30163 K. For a H 2 rate oc exp(-y/rki n ), the exponential fac- 
tor for H* 2 is thus y* = max(0, y - 30163K). This approximation 
however overestimates the rate coefficients, if y is too large and 
we switch off reactions with y > 20000 K. 

The population of H^ is determined by collisions with H and 



H 2 (Tielens & Hollenbach 1985), spontaneous decay (2 x 10 



Londonll 19781)" and the pumping by FUV photons, which is 



assumed to be ten times the photodissociation rate. 
A.4. Molecular/atomic excitation 

The population n, (in cirT 3 ) of a level i of a molecular or atomic 
species is determined by the rate equation 

^ = X njPtj - m J] Pij + Ft ~ A = 0, (A. 10) 

with the rate coefficients for the transition between level i and j, 



Aij + BijiJ^ + dj (/•:, > /•:,) 



Bij(Jij) + Cij 



(Ei < Ej) 



(A.ll) 



with the level energy the collisional rate coefficients Cy and 
the Einstein A and B-coefficients Ay and B/j, respectively. The 
rate of chemical formation of a molecule forming into the state i 
and destroyed from that s tate is given by F{ and £),-, respectively 
( Staub er & Brudereii 2009). In the following, we assume steady- 
state conditions and set dni/dt = 0. To calculate the ambient 
radiation field 



(Jij) = ^ff <Pij(y)I v (Q.)dvdn , 



(A. 12) 



the normalized line profile function 0y(v) an d the intensity de- 
pending on the frequency and direction is needed in every cell 
of the model. The intensity along a ray and for one frequency is 
calculated from the radiative transfer equation 

dl v 

— = -a v I v + j v , (A. 13) 

ds 

with the absorption and emission coefficient a Y and j Y , respec- 
tively and the distance ds that the ray propagates. Both coef- 
ficients contain contributions from both line and dust absorp- 
tion and emission. Since the level population enters the ab- 
sorption and emission coefficients, Eqs. I A. 101 IA.12I and IA. 131 
form a coupled problem, whic h is computationally very demand- 
ing to solve (see for exampl e iHogerheiide & van der Sl|2000; 
iBrinch & Hoger heiide 201QI). 

In this work, we employ an approximation similar t o the es- 
cape probability method (e.g. Ivan der Tak et a fl 120071) . In this 
approximation, the physical conditions and thus the source func- 
tion S v = j v /a v is assumed to be constant in the modeled region. 
The ambient radiation field (J/j) can then be expressed by an an- 
alytical expression, which greatly facilitates the calculation. For 
our models with physical conditions strongly varying with posi- 
tion, the approximation of a single zone with constant physical 
conditions is not applicable. The solution of Eq. lA.13l then reads 



S v(t' v )< 

Jo 



dj' 



(A. 14) 



using the definition of the optical depth, dr v _ a v ds and th e 
background intensity 7 v> b g . Similar to lPoelman & Spaansl (2005), 
we solve for the excitation at every position, although we ap- 
proximate the local radiation field by keeping the source func- 
tion constant along a ray. The time-consuming integration in Eq. 
IA. 141 reduces in this way to a simple calculation of the opacity 
along a ray from the current position to the edge of the mod- 
eled region. This approximation allows us to solve for the exci- 
tation of complex problems (many lines, high optical depth) in 
a short time. For example the excitation of water in the models 
presented in this work can be solved in a few minutes using a 
standard PC. 

To include pumpin g by dust and the vel ocity structure, we 
follow the approach of Taka hashi et al.l (Il983l) . For the dust and 
line opacity t Vi l and t Vi d, along a ray, we calculate 



j_ rr t v ,d + r VtL e 

4?r J J t v>d + T V ,L 



(Tv,L+Tv,D) 



t>ij(v)dvdD. 



-(t,,l+t v . d ). 



>ij(v)dvdCl . 



(A.15) 



(A. 16) 



(A.17) 



The level population is then obtained from Eq. lA.lOl with 

_(A ij e ij + B ij (J' ij ) + C ij (Ei > Ej) 

U ~ I BijiJlj) + Cy {E t < Ej) , 

and (J'..) = (ey - r/ij)B v (T D ) + rjjjI Vy \, g . The local dust temperature 
is given by To and the background radiation (CMB) is given by 
Iv,bg- For the calculation of ey and jyy, either 6 or 26 directions 
and of order 40 frequency bins per line are used. 

Once the level population is determined, the molecular cool- 
ing rate can be calculated and synthetic maps can be derived 
by solving the radiative transfer equation (Eq. IA.14I ). Our ray- 
tracer produces images in fits format. To properly resolve the 
inner structure (hot core or inner disk), we shoot several rays per 
pixel through that region and average over the intensity. Along 
the ray, the integration steps within one computational cell are 
refined, according to the velocity structure in order to properly 
resolve narrow lines. 
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Benchmark test 

To benchmark th e excitation calculation, we calculate the prob- 
lem proposed by Ivan Zadelhoff et a! (2002), which consists of 
the calculation of the HCO + excitation in a spherical, infalling 
class cloud core. The physical conditions are given in Fig. IA.5l 
The abundance of HCO + is fixed to either 10~ 8 or 10~ 9 relative 
to the H2 density. The difficulty of this benchmark problem lies 
in the high optical depth reached in the lines and that collision 
partner densities are below the critical density. The level popula- 
tion is thus far from the local thermal equilibrium and radiative 
excitation is important. 




2 ■ 10 10 5 ■ 10 16 1 ■ 10 17 2 ■ 10 17 
Radius (cm) 

Fig.A.5. Physical conditions of the cloud core used to bench- 
mark the excitation calculation. 



In Fig. IA.61 we compare the solution obtained with 
this method with the solutio n calculated by RATRAN 
(Hogerheiide & van der Tak 2000). That code is widely used and 
has been well-tested against other codes. For an abundance of 
10~ 9 , the normalized level population agrees to within 30% with 
the solution provided by RATRAN. An exception is the J — 3 
level in the inner part of the cloud. The deviations are larger for 
an abundance of 10~ 8 , but still mostly within 30%. The intensity 
is obtained with the raytracer implemented in our code and con- 
volved to the telescope beam. In Fig. IA.61 the / = 1 — > and 
J = 4 — > 3 lines are shown for a beam of 29" and 14", respec- 
tively. We also indicate the intensities obtained assuming local 
thermal equilibrium (LTE) and "thin" conditions, setting the am- 
bient radiation field to 0. As the level populations, the line fluxes 
mostly agree to within about 30% of the intensity obtained with 
RATRAN. We note that the intensities obtained with RATRAN 
are calculated with their raytra cer (SKY) and convolved using 
the MIRIAD package dOill2005l) . Thus, this benchmark problem 
also tests our raytracer and convolution routine. 

The agreement fo und in with this bench mark problem is sim- 
ilar to the findings o f IWoitke et alJ <|2009b), who calculated the 
water emission from a Herbig Ae disk calculated with a simi- 
lar method and compared them to RATRAN. We consider the 
agreement as reasonable and sufficient for our application, since 
the uncertainties entering through the chemical network calcu- 
lation are much larger. It also shows that this method is a much 
better approximation than e.g. assuming LTE. 



A.5. Thermal balance 

The gas temperature is obtained from the equilibrium between 
the heating and cooling rates 

i i 

where e is the internal energy of the gas and T, and A, are dif- 
ferent heating and cooling rates (for example, line cooling or 
photoelectric heating on dust grains). This non-linear equation is 
solved by a secant method starting from the dust temperature. If 
the convergence is insufficient, the process is restarted with a dif- 
ferent value. Convergence is reached if the heating and cooling 
rate agree to within a predefined threshold in all cells (typically 
1%). We implement th e following heating and cooling rates (see 
iBruderer et all2009al) : 

1. ) Photoelectric heating 

Photoelectric heating o n large graphite and s i licate grains is 
implemented following iKamp & van Zadelhoftl d200lT) . The re- 
quired grain absorption cross-section, averaged over the FUV 
band, is taken to be consistent with the calculation of the UV 
field and the dust temperature. As for the chemistry, the rates 
involving dust grains are scaled to the local gas/dust ratio. 
Photoe lectric heating on small dust grains is implemented fol- 
lowing Bakes & Tielensl(ll994l) . We also implement their recom- 
bination cooling rate. 

2. ) Gas-grain heating or cooling 

Depending on the difference between gas and dust temperature, 
gas-grain collisions can either heat or cool the gas. In regions 
of very high density, gas-grain colli sions couple the gas tem per- 
ature to the dust temperature (e.g. iDotv & Neufeldlll997l). We 
imple ment the heating and cooling rates following I Young et alJ 
(I2004I) . 

3. ) H 2 heating 

Molecular hydrogen can contribute to the heating and cooling 
in different ways: (i) Line cooling, (ii) Heating through colli - 
sional deexcitation of FUV pumped H2. (Hi) Formation heat- 
ing, (iv) Photodiss ociation heat i ng. F or (i) and (ii), we use 
the analytica l fit of Rollig et al. (2006) to the exact multilevel 
treatment by Sternberg & Dalgarno] d 1995b . For (Hi), we follow 
IKamp & van Zad elhoff (2001)" and assume that one third of the 
binding energy of ~ 4.5 eV is released to the gas for heating. 
Photodissociation (i v) carries away abou t 0.4 eV in the form of 
heat to the gas (e.g. iJonkheid etal]|2004l) . We note that the H 2 
dissociation rate is taken to be consistent with the chemistry. 

4. ) Cosmic-ray heating 

Follow ing ICravens & Dalgarnold 19781) andfGlassgol d & Langerl 
(Il973l) . a cosmic-ray ionization of H 2 (H) releases 8 eV (3.5 eV) 
to the gas. 

5. ) X-ray heating 

We im plement the X-ray heating rates following Dalgarno et al. 
(1999) for a H 2 , H, and He mixture. The energy deposition rate 
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Fig. A.6. Result of the benchmark of the excitation calculation. The upper panel shows the results for a HCO + abundance of 10~ 9 
relative to H2, the lower panel for an abundance of 10 relative to H2. The plots on the left display the normalized level population of 
the first eight levels, the plots on the right show beam-averaged intensities. The black lines give results obtained using the RATRAN 
code, and the red dots/lines show results derived by this work. The gray shaded area indicates a 30% range to the RATRAN solution. 



per particle Hx is calcu lated using the cross-sections provided in 
iBruderer et ail d2009bl) . 



8.) Carbon ionization 



6.) Line cooling (molecular and atomic fine-structure lines) 



Cooling rates of molecular lines and atomic fine-structure lines 
are derived from the level population calculation (Sect. IA.4I) , 
We account for O, C, C + , CO, 13 CO, H 2 0, and OH. The abun- 
dance of 13 CO is obtained from scal ing the CO abundance by 
12 C/ 13 C = 70 dWilson & Roodll 19941). Molecula r data are taken 
from the LAMDA database ( Schoi er et al.ll2005h and include the 
collisional excitation with H, H2, and electrons for O, C, and C + , 
respectively. For the molecular species, only excitation with H2 
is available. 



7.) Hydrogen Lya- and [O I] 6300 A cooling 



Hydrogen Lya cooling and cooling by neutral oxygen [O I] by 
means of the metastable 'D- 3 P line at 6300 A is included fol- 
lowing ISternberg & Dalgarnoldl989l) . 



Carbo n ionization delivers about 1 eV to the gas dJonkheid et alj 
2004). The ionization rate is assumed to be consistent with the 
chemical network including the effects of self-shielding. 

A.5.1. Benchmark test 

To test the implementation of the thermal balance calculation, 
we run the benchm ark problem of the PDR comparison study of 
Rolli g et al] (I20071) . which consists of four slab models with a 
density of 10 3 and 10 5 ' 5 crrT 3 and an incident UV radiation fiel d 
of x = 10 or x = 10 5 (in units of Draine ISRF, lDraindll978b . 
To eliminate the differences caused by the adopted chemical net- 
work, the benchmark study is run with the same simple chemical 
network consisting of only the 3 1 species used in th e PDR com- 
pariso n study. The parameters given in Table 5 of Rollig et al. 
d2007l) are implemented. 

The gas temperature calculated with our code is shown in 
Fig. IA.7l together with the results from other PDR models. The 
gas temperature derived with our code is in reasonable agree- 
ment with that found by other PDR codes. Differences between 
the different codes are however considerable, particularly for 
the high-density (10 55 cm" 3 ) model with high UV intensity 
(X = 10 5 ). This model, is however closest to the conditions in 
disk atmospheres and outflow walls. 
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a.) Gas- Temperature (n = 10' f cm , \ = 10) b.) Gas-Temperature (n = 10' f cm ,! , \ = 10 :) ) 
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Fig. A.7. Result of benchmark test for the calculation of the gas temperature. The results of the codes participating in the study by 
Rol hgetal] (120071) are shown in gray, and the gas temperature obtained with our code is given in red. 



Appendix B: Analytical estimates of fluxes 

B.I. Analytical estimate of the [CWJflux 

To analytically estimate the [C II] flux that might originate from 
a remnant envelope or the foreground, we use 



A hv 

F = dn-I=-z—A ul N c *x u , 
d z An 



(B.l) 



for the flux F, the integrated intensity /, and the solid angle dQ.. 
The solid angle of the emitting area is calculated from the pro- 
jected area A and the distance d to the source. The integrated 
intensity / is obtained from the line frequency v, the Einstein-A 
coefficient A„i, the column density of C + Nc+, and the normal- 
ized level in the upper level x u . Assuming the level population 
in the LTE, which is the case above the critical density of a few 
times 10 3 cm -3 , we reach a maximum x u = 2/3, given by the 
statistical weights of the C + levels. These equations assume op- 
tically thin emission. For a line width of 1 km s , r = 1 at the 
line center is reached for a column density of a few times 10 18 



B.2. Rotation diagram of CO 

For the rotation diagram of CO, we consider the optically thin 
flux of CO, which is assumed to be thermalized, 



-EJkT 



da- 1 



ui 



dn^A ul N(CO)^— 
4n Q(T) 



(B.2) 



using the same notation as in Appendix lB . 1 1 and the CO column 
density N{CO), the upper level energy E u , and the partition func- 
tion Q{T). Rearranging yields 



e Y = 



4jtF„ 



= dON(CO)- 



-EJkT N 



Auihv u ig u ' Q{T) 
with the column density per level N u = N{CO) ( , T , 



(B.3) 



Thus, 



Y = In 



4ttF u 



A u ihv„ig u 



InU^))-^ 
I Q(T) ) kT 



(B.4) 



We note that the rotation diagram is often given by assuming 
that F„i = dClbeamluh for the solid angle of the telescope beam 
dQ., thus the molecules are equally distributed over the beam. 
This is however inappropriate for the CO ladder probing a wide 
range of temperatures and subsequently different emitting re- 
gions (Section l3~4b . 

How do opacity effects alter the rotation diagram? Assuming 
a square line profile of (velocity) width Av for simplicity and 
neglecting the background intensity, we obtain for the integrated 
intensity using Eq. IA.141 >. 



I ul = Av—B(T ex , ul )(l-e- T "<) 
c 



(B.5) 



For t u i <K 1, /„/ ~ Av-^ i Z?(r eXi „;) T «; = jvdl, recovering the ex- 
pression used in Eq. IB.2l For r„/ > 1, /„/ ~ Av^-B{T ex u i), which 
is always larger than the optically thin expression. Thus, opti- 
cal depth effects shift points above the optically thin lines in the 
rotation diagram. 



